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1  INTRODUCTION 


l.I  Background 

This  project  was  a  team  effort,  originally  scheduled  two  years,  from  July  2001  to  June  2003,  then 
extended  till  12/31/2003,  led  by  CFD  Research  Corporation  (CFDRC),  with  participation  of 
University  of  California  Santa  Cruz  (UCSC),  Cornell  University  (CU),  and  University  of 
Delaware  (UD). 

The  objective  of  this  project  was  to  develop,  demonstrate,  and  validate  Design  Tools  for  Opto- 
Electronic  Integrated  Circuits  and  Interconnects.  The  new  tools  include  several  innovative  ideas 
in  the  following  Technical  Areas: 

•  EDA  Interfaces  and  Front  End  for  Automation, 

•  Efficient  Extraction  of  Parasitics  in  Electrical  Interconnects, 

•  New  Methods  of  Electrical  Interconnects  Modeling  for  High  Frequencies, 

•  New  Methods  of  Optical  Interconnects  Modeling, 

•  Model-Order  Reduction  and  Behavioral  Models  for  Electronic-Photonic  Devices  and 
Circuits, 

•  Integration  of  the  New  Methods  and  Software  into  EDA  Tools  for  Opto-EIectronic  (0/E) 
Design  Flow  and  Automated  Synthesis. 

Major  activities  in  this  program  were  organized  to  take  advantage  of  the  expertise  of  all  the  team 
members  and  their  respective  Pi’s: 

•  CFDRC  (Dr.  Marek  Turowski,  Dr.  Andrzej  Przekwas,  Dr.  Robert  Tramel,  Dr.  H.Q 
Yang):  Computational  Microelectronics,  Semiconductor  Physics,  Electromagnetics, 
Photonics,  0/E  Interconnects,  Computational  Physics,  Appl.  Math.;  Numerical  Methods 
and  Algorithms,  Multi-Disciplinary  CAD,  TCAD,  project  coordination,  software 
development,  validation,  demonstration  and  applications. 

•  University  of  California  Santa  Cruz  (Dr.  Steve  Kang):  VLSI  CAD;  High-performance 
VLSI  Systems;  Modeling  and  simulation  of  high-speed  electronic  and  optoelectronic 
integrated  circuits  and  systems;  Reduced/behavioral  models  of  0/E  devices  and  systems. 

•  Cornell  University  (Dr.  Edwin  Kan):  Nanoscale  electronic  devices  and  interconnects; 
Technology  CAD  and  intelligence  abstraction  for  VLSI  systems;  Computational 
quantum-electronics  for  quantum-well  photonic  devices. 

•  University  of  Delaware  (Dr.  Fouad  Kiamilev):  Design  and  Test  of  high-performance 
Digital  and  Analog  CMOS  and  BiCMOS  VLSI  Circuits  for:  Photonics,  Optoelectronic 
Integrated  Circuits  (OEIC),  and  Optoelectronic  VLSI;  RF  IC  Design,  Fabrication  and 
Measurements  for  Verification  and  Validation  the  newly  developed  software. 

This  report  documents  the  work  performed  during  the  30  months  of  the  project. 
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1.2  Project  Motivation 


Decreasing  feature  sizes  of  very  deep  submicron  ULSI  (Ultra  Large  Scale  Integration)  as  well  as 
increased  operating  frequencies  of  mixed-signal  (digital/analog/RF)  and  mixed-technology 
(electronics,  photonics,  etc.)  systems  make  it  necessary  to  precisely  predict  parasitic  as  well  as 
intended  electromagnetic  (EM)  and  optical  effects  of  interconnects.  Equally  important  becomes 
the  accurate  simulation  of  the  interaction  between  active  devices  (MOSFETs,  VCSELs,  etc.)  and 
interconnects. 

High-frequency  and  high-bandwidth  circuits  can  only  be  achieved  through  nano-scale  devices 
due  to  the  upper  limit  of  information  propagation  speed.  Above  30GHz,  the  wavelength  of  EM 
waves  becomes  1cm,  and  the  distinction  between  processing  unit  (on  the  order  of  2cm  chip)  and 
transmission  units  (antenna  or  optical  fiber  connection)  becomes  blurred.  To  support  new 
system  design  in  this  domain,  not  only  the  CAD  system  needs  to  deal  with  the  ever-increasing 
complexity  and  component  interconnect  coupling,  but  also  new  nano-scale  effects  needs  to  be 
included  in  a  more  rigorous  and  universal  way.  It  is  then  possible  to  integrate  the  nonlinear 
characteristics  into  the  system  design  validation  and  optimization. 

Limitations  of  the  methods  and  tools  existing  before  NeoCAD  Program: 

-  Few  EM  solvers  applicable  for  interconnects  (no  automatic  model  building,  limited  model 
area). 

-  No  existing  inductance  extraction  tools  (identifying  low-impedance  return  path  requires 
analysis  of  the  global  circuit,  which  was  prohibitively  compute-intensive  with  the  existing 
field  solvers). 

-  Equivalent-circuit  models  of  ULSI  interconnects,  used  with  SPICE,  become  computationally 
inefficient  and  often  prohibitive  because  of  their  size  and  lack  of  circuit  matrix  sparsity  due 
to  distributed  nature  of  inductive  and  capacitive  coupling  (Figure  1) 

-  No  available  optoelectronic  (0/E)  electric-optical  device  and  system  simulators. 

-  No  coupled  Device+EM  modeling  capability,  e.g.  for  RF  SiGe  HBT,  BiCMOS,  HEMT,  etc. 
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Figure  1.  SPICE-type  equivalent-circuit  models  of  ULSI  ICs,  interconnects,  and  mixed-signal  RF 
circuits  become  computationally  inefficient  and  often  prohibitive  for  higher  frequencies. 
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Need  for  relevant  TCAD  and  EDA  tools  for  opto-electronic  (0/E)  circuits  and  systems 

The  increasing  design  complexity  of  the  current  and  future  ICs  on  each  design  level  (component, 
module,  package,  system)  is  enhanced  even  more  by  introduction  of  optical  elements  on  chip  and 
printed  circuit  board.  This  calls  for  new  modeling  and  simulation  approaches  allowing  a  fast  and 
accurate  interconnect  modeling  where  a  comprehensive  treatment  of  all  relevant  physical  effects 
has  to  be  considered.  Simply  speaking,  such  comprehensive  simulation  and  design  tools  for  0/E 
circuits  and  systems  do  not  exist  yet.  Some  of  the  reasons  are  as  follows: 

■  the  comprehensive  analysis  of  0/E  circuits  require  coupled  simulations  in  mixed  domains, 
like  electrical,  mechanical,  thermal,  and  of  course  full-vector  electromagnetic; 

■  three-dimensional  (3D)  and  spatially-dependent  behavior  of  optical  systems  (Fig.  16,  12,  and 
15)  can  not  be  directly  captured  by  standard  EDA  tools  based  mostly  on  2D-layout  approach; 


Figure  2.  Example  of  a  hybrid  3D  opto-electronic  design  with  a  communication  waveguide. 

■  the  wave  nature  of  optics  and  high-frequency  electronics  requires  wave-based  and  vector- 
based  description  and  modeling,  which  makes  SPICE-type  equivalent-circuit  modeling  very 
difficult  or  sometimes  impossible; 

■  no  comprehensive,  reliable  behavioral  or  equivalent-circuit  models  of  photonic  devices  and 
next-generation  electronic  compound  heterojunction  devices; 

■  no  available  electric-optical  device  simulators  for  new  optoelectronic  devices  (see  below). 

Why  the  existing  TCAD  tools  could  not  efficiently  handle  0/E  devices... 

Before  NeoCAD,  there  was  no  commercial  TCAD  tools  on  the  market  which  would  handle  all 
the  complexity  of  the  photonic  device  operation,  such  as  for  example  VCSEL.  No  existing 
device  simulator  contained  fully  coupled  Device  +  Electromagnetic  modeling  capability  (for 
accurate  light  propagation  modeling),  enhanced  with  satisfactory  model  of  quantum  well  physics 
as  well  as  thermal  phenomena  influencing  both  electrical  and  optical  behavior.  The  only  existing 
tools  on  the  market  had  only  partial  components  for  optical  devices  simulations,  like  for  example 
Crosslight  software  from  Canada,  or  ISE  TCAD  tools  from  Switzerland. 

1.3  Project  Objectives  and  Scope 

The  overall  objective  of  the  OE-IDEAS  project  was: 

*X*  Design  Automation  Tools  for  Opto-Electronic  Integrated  Circuits  and  Interconnects 
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Our  intension  was  to  develop,  implement,  and  demonstrate  several  innovative  ideas  in  the 
following  Teehnical  Areas: 

■  EDA  Interfaees  and  Front  End  for  Automation 

■  Effieient  Extraetion  of  Parasities  in  Eleetrieal  Intereonnects 

■  New  Methods  of  Eleetrieal  Interconnects  Modeling  for  High  Frequencies 

■  New  Methods  of  Optical  Interconnects  Modeling 

■  Model-Order  Reduction  and  Behavioral  Models  for  Electronic-Photonic  Devices  and 
Circuits 

■  Integration  of  the  New  Methods  and  Software  into  EDA  Tools  for  Opto-Electronic 
Design  Flow  and  Automated  Synthesis 


The  scope  of  this  project  in  the  overall  Design  Flow  Methodology  for  mixed-signal  electronic- 
photonic  circuits  is  illustrated  in  the  figure  below. 
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Figure  3.  Extended  Design  Flow  of  electrical-optical  circuits  and  printed  circuit  boards, 

with  indicated  scope  of  this  project. 


As  microelectronics  technology  continues  to  advance,  the  associated  electrical  interconnection 
technology  is  not  likely  to  keep  pace,  due  to  many  parasitic  effects  appearing  in  metallic 
interconnections.  Hence,  the  optical  interconnections  and  photonic-electronic  chips  seem  to  be 
the  most  viable  technology  for  future  multi-GHz  mixed-signal  electronics  [Li  2000].  Therefore, 
this  projects  focused  on  the  following  crucial  issues  in  this  domain: 

1.  Efficient  and  fast  techniques  of  extraction  of  parasities  as  well  as  revolutionary  new 
methods  of  modeling  parasities  in  electrical  interconnects,  where  high-frequency 
electromagnetic  effects  will  not  allow  for  using  SPICE-type  models  anymore 
(cf.  Figure  1); 
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2.  New,  comprehensive  CAD  tools  for  integrated  opto-electronic  (0/E)  circuits  and 
systems,  which  would  allow  designers  to  analyze  and  optimize  the  full-path  signal 
characteristics,  at  all  the  stages  of  the  optically  communicating  circuits  (Figure  4). 
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Figure  4.  Examples  of  photonic  VLSI  integrated  circuits  and  interconnections  for  future  high- 
bandwidth  communication  systems,  developed  under  DARPA  FSOI  and  VLSI  Photonics  programs. 


Such  a  comprehensive  0/E  CAD  tool  will  allow  to  obtain  the  entire  signal  characteristics  of: 

■  “Electrons  to  Photons”  Path:  from  electrical  signal  drivers  through  all  the  complex 
interconnects  and  bonds  of  multi-chip  setup  up  to  optical  wave  characteristics  at  the 
output; 

■  Optical  Path:  through  lenses,  waveguides,  fibers,  resonators,  etc. 

■  “Photons  to  Electrons”  Path:  from  optical  flux  input  through  all  the  complex 
interconnects  and  bonds  down  to  electrical  analog-to-digital  (A/D)  converters  and 
amplifiers. 


1.4  Summary  of  Accomplishments 

All  of  the  accomplishments  aimed  for  fulfillment  of  the  main  objective  of  this  project  and  its 
components,  as  described  in  the  previous  section.  Following  subsections  list  the  most  important 
accomplishments. 

•  First  in  the  U.S.  commercial  software  for  high-fidelity  multi-mode  simulation  of 
Lasers  (VCSELs,  edge-emit.)  and  Photodetectors,  O’SEMI  from  CFDRC: 

-  new  Optical  Gain  Model,  with  quantum-physics  and  thermal  phenomena; 

-  implemented  Effective  Index  Method  for  Edge-Emitting  Laser  optics; 

-  new  Tunneling  Current  Model,  based  on  quantum  model  of  heterojunction; 

-  currently  installed  and  used  by  AFRL/WL,  ARL  (Adelphi,  MD),  NRL,  UCSB, 

UC  Santa  Cruz,  U.  Florida,  U.  Arizona,  U.  Waterloo,  Vanderbilt  Univ.,  and  others; 

-  inquiries  from:  Sandia  National  Labs,  LLNL,  McMaster  University,  et  al. 

•  New  capabilities  in  CFDRC’s  Full-Wave  Electromagnetic  Simulators 
for  RF  and  Optical  Interconnects: 

-  Signal  monitoring  and  postprocessing:  V(t),  I(t),  Z(w),  S-param(w); 

-  Electromagnetic  filament  model  implemented  and  initially  tested; 

-  CFDRC  full-wave  EM  solvers  offered  to  NeoCAD  Community,  and 
used  by  Univ.  of  Delaware,  U.  of  Washington,  and  U.  of  Cincinnati. 

.  CFDRC’s  Full-Wave  3D  Simulator  EMAG  coupled  with  SPICE; 
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•  Demonstrated  Photonics  Applications  and  provided  Tutorials  for  EMAG  users: 

1)  Optical  waveguide  coupling; 

2)  Optical  ring  resonator; 

3)  Photonic  crystal  waveguide; 

4)  Photonic  bandgap  2D  matrix. 

•  Simulated  several  NeoCAD  benchmark  structures: 

1)  CPI  waveguide  from  HRL-Mayo, 

2)  20G  transmission  line  from  NCSU, 

3)  MIT-CLineOl  coplanar  waveguide, 

4)  NCSU-2BA  Coupled  Spiral-Inductor  Transformer, 

5)  100  GHz  100  GHz  InP  circuit  from  Mayo  (see  illustration  below). 


Figure  5.  Parasitic  Extraction  of  100  GHz  InP  circuit,  a  NeoCAD  benchmark  from  Mayo. 


CFD-Micromesh  -  a  unique,  automated  3D  mesher  from  IC  layouts  -  enhanced,  offered 
to  NeoCAD  Community;  new  tutorials  provided. 

Tasks  led  by  University  of  Delaware: 

-  Designed,  fabricated,  tested  and  measured  an  optical  transceiver  chip,  and  built  a 
3D  model  (with  CFD-Micromesh). 

-  Characterized  and  built  a  3D  model  (with  CFD-Micromesh)  of  the  Mayo  lOG  system 
multi-chip  module  (MCM). 

-  Formulated  the  transmission-line  model  and  the  7i-model  of  via-stub  resonance  on 

system  PCB’s;  with  CFDRC  simulated  the  model  performance  . 
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-  Designed,  fabricated,  and  tested  Dual-Rate  Current  Mode  Logic  electrical  10 ’s. 

•  Tasks  led  by  Cornell  University: 

-  Established  the  quantum  transport  model  hierarchy  in  density  gradient,  quantum 
hydrodynamic,  full  Wigner,  and  NEGF  on  large  barrier  tunneling  problems.  The 
benchmark  problems  are  RTD  and  MOSFET  gate  tunneling. 

-  Established  large-signal  transient  analysis  from  S-parameters  on  nonlinear  circuits. 
The  benchmark  problems  include  high-power/high  frequency  transistors,  nonlinear 
transmission  lines,  and  buffer-insertion  interconnect. 

-  Proposed  and  analyzed  pulse-based  interconnect  architecture  and  circuit 
implementation,  based  on  the  large-signal  transient  analysis. 

-  Designed  and  verified  the  use  of  patterned  permalloy  to  control  the  on-chip  magnetic 
fields  for  accurate  partial  inductance  extraction. 

-  Designed  and  verified  the  use  of  undergrids  to  control  CMOS  interconnect 
dispersion  on  transmission  lines  and  co-planar  waveguides. 

•  Tasks  led  by  UC  Santa  Cruz: 

Developed,  implemented,  and  tested  compact  models  for  OE  devices  for  mixed-level 
mixed- technology  simulation,  and  demonstrated  in  VHDL-AMS  (Figure  6). 

Modeling  VCSEL  and  transmitter: 

-  Developed  a  behavioral  VHDL-AMS  model  of  VCSEL; 

-  Implemented  the  model  into  SMASHTM  from  Dolphin  Integration; 

-  Tested  the  VCSEL  model  by  time  domain  simulation  of  different  devices; 

-  Simulated  a  simple  transmitter  consisting  of  CMOS  driving  circuit  and  VCSEL. 
Modeling  Metal-Semiconductor-Metal  photodetector: 

-  Developed  a  behavioral  VHDL-AMS  model  of  MSM  photodetector; 

-  Implemented  the  model  into  hAMster™  provided  by  AnSoft; 

-  Validated  the  model  by  comparing  simulation  results  with  experimental  data. 

Modeling  Raman  Amplifier: 

-  Developed  a  behavioral  VHDL-AMS  model  of  Raman  Amplifier; 

-  Implemented  the  model  into  hAMster  from  AnSoft; 

-  Validated  the  model  by  comparing  simulation  results  with  experimental  data. 

Modeling  Erbium  Doped  Fiber  Amplifier  (EPF A): 

-  Developed  a  behavioral  VHDL-AMS  model  of  EDFA; 

-  Implemented  the  model  into  hAMster  from  AnSoft; 

-  Validated  the  model  by  comparing  with  the  experimental  data. 

Automation  of  0/E  behavioral  model  parameters  extraction: 

-  Finished  the  linear  part  of  VCSEL  behavioral  model  parameter  extraction  using  Python; 

-  Implemented  a  nonlinear  behavior  model  parameter  extraction  algorithms  for  Python. 
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High-Fidelity  CFDRC  Models  Behavioral  VHDL-AMS  Models 


Figure  6.  Extraction  of  OE  behavioral  models  for  mixed-level  system  simulation. 

2  OPTO-ELECTRONIC  DEVICE  MODELS  CONSOLIDATION  FOR  BEHAVIOR 
MODEL  EXTRACTION 


2.1  Enhancement  of  High-Fidelity  Electro-Optical  VCSEL  Model 

During  this  project  CFDRC  was  working  on  enhancement  of  the  high-fidelity  model  of 
semiconductor  devices  to  take  into  account  the  complex  phenomena  inside  Vertical-Cavity 
Surface-Emitting  Lasers  (VCSELs).  This  will  allow  for  more  realistic  modeling  and  simulation 
of  VCSELs  and  other  photonic  devices,  also  for  better  quality  compact  model  extraction. 

The  new,  enhanced  electro-optical  model  includes  the  following  equations: 


Electric  Potential  (Poisson ’s)  Equation 
V-{-eV(^)  =  q[p-n  +  Nt)-N-A) 


Carrier  Continuity  Equations 


q--VJ„=q{a-R) 

ot 

q^^+V-Jp=q{G-R) 

where  the  current  densities  are: 


J„  =  qDyn  +  qn 


r 

(  e'] 

3  1 

1  ^  J 

+DyT„ 

High-Fidelity  VCSEL  Model 
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Jp=-qDpVp  +  qp\ 

Carrier  Energy  Equations 


r 

+ 

1 

dnw, 

dt 

Spwp 


^  +  v-s„-v 


-0  +  - 


^  J 


^wn 


+  V-S„-V 
dt  P 


\ 


V  ^  J 


W„-Wt  ,  , 

■Jn+P^ - -^wJG-R) 


'"wp 


where  energy  fluxes  are 

Sn  ^-kn^Tn--  i^n  +  kTn  )j „ 
q 

Sp  =  —kpVTp  A  {w p  +  kTp ^ p 

q 

Photon  Number  (Rate)  Equation 

^Sm  ~  2p  ^  J  ^sp^^  ~  ® 


Optical  Wave  Equation 
v^w+k^[s-jq^}v  =  o 


Example  results  of  Quantum  Well  (QW)  modeling  are  shown  in  the  figures  below.  Such  QW 
structures  are  typical  for  photonic  devices,  including  VCSELs.  The  n  and  p  regions  at  the  ends  of 
the  structure  in  Figure  7  represent  lumped  models  of  the  Distributed  Bragg  Reflectors  (DBR). 


n  -  AIo.625G3o.375AS 


AI0.25G30.75AS 


hUi 


P  ■  Alo.625Gao.375AS 


spacer 

Alo.8Gao.4As 


TTT 

GaAs 


Bias  2V 


Figure  7.  High-fidelity  multiple  Quantum  Well  (QW)  model. 
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Figure  8.  Energy  Band  Diagram  of  the  Multiple-QW  structure. 


Figure  9.  Carrier  Density  near  the  Multiple-QW  structure. 


The  results  shown  above  include  the  following  new  OE-IDEAS  developments: 

•  Heterostructure  Model  (Energy  Band  Gap,  Thermionic  Current) 

•  Quantum  Well  Model  (Tunneling  Current) 

•  Full  Newton  Iterative  Algorithm 

2.2  Coupled  Multiphysics  Modeling  of  Semiconductor  Lasers 

Comprehensive  models  of  semiconductor  lasers  are  required  to  predict  realistic  behavior  of  various 
laser  devices  for  the  spatially  nonuniform  gain  that  results  due  to  current  crowding.  Nonuniform 
gain  has  visible  effect  on  laser  dynamics  and  parameters,  like  threshold  gain.  Consistent  solution  of 
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coupled  electrical,  thermal,  gain  and  optical  problems  was  achieved  by  CFD-ACE+  integrated 
solver.  Integration  of  multiphysics  into  a  single  computational  environment  allowed  for  high- 
efficiency,  high-fidelity  modeling  of  modem  semiconductor  lasers,  including  vertical-cavity 
surface-emitting  lasers  (VCSELs)  and  edge-emitting  lasers  (EELs). 

This  work  allowed  for  the  first  time  a  full  self-consistent  coupling  of  the  most-advanced-physics 
models  of  spatially  dependent  current  flow,  temperature  effects,  gain,  and  optical  modes 
development  and  competition.  Three  coupled  modules  of  CFD-ACE+:  Semiconductor  Device, 
Thermal,  and  Optics  are  used  in  the  present  study.  The  Semiconductor  Device  module  solves 
electric  model  based  on  drift-diffusion  (DD)  equations  for  carrier  densities,  and  energy  balance 
(EB)  equations  for  carrier  temperatures.  It  has  many  advanced  capabilities  required  by  VCSEL 
simulations,  such  as  quantum  well  (QW),  heterostmcture  interface,  etc.  Figure  10  shows  an 
example  of  simulation  results  of  current  flow  and  electron  density  distribution  for  an  intracavity 
VCSEL  problem. 

Optical  models  include  Weighted  Index  Method  (scalar)  and  Method  of  Lines  (scalar  and 
vectorial).  Both  edge-emitting  lasers  and  VCSELs  can  be  efficiently  modeled.  The  code  allows  to 
study  transient  behavior  of  lasers  with  time-dependent  and  spatially  inhomogeneous  gain.  Multiple 
quantum  wells  and  multiple  transverse  modes  are  included,  as  illustrated  in  Figure  11  through 
Figure  14.  For  each  time  step,  a  calculation  loop  over  electrical,  thermal,  gain,  and  optical  solvers 
is  performed. 

Solution  for  VCSEL  transverse  and  axial  field  distribution  in  the  modes  as  well  as  their  lasing 
frequencies  consists  of  solution  of  an  eigenvalue  problem  for  Helmholtz  equation  in  cylindrical 
coordinates,  and  boundary  value  problem  set  by  the  boundary  conditions  at  VCSEL  layer  and 
radial  region  interfaces.  For  each  mode  photon  rate  equation  is  solved  including  optical  losses,  and 
spontaneous  and  stimulated  emission  terms  taking  into  account  spatially  and  time-dependent 
material  gain.  Corresponding  stimulated  and  spontaneous  emission  terms  are  responsible  for 
carrier  relaxation  in  the  electrical  model. 

The  simulation  results  were  validated  through  comparison  with  experimental  data  obtained  from 
U.S.  Army  Research  Laboratory  (ARE,  Adelphi,  MD).  This  kind  of  modeling  is  crucial  for 
comprehensive  analysis  design,  and  optimization  of  new  architectures  of  semiconductor  lasers,  as 
well  as  for  reliable  extraction  of  compact  model  parameters  for  system-level  simulations. 
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Figure  10.  Current  flow  and  electron  density 
distributions  inside  an  intra-cavity  VCSEL 
device 
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Figure  12.  Steady-state  radial  distribution  of 
intensity  in  different  transverse  modes  for 
VCSEL  with  6  micrometer  oxide  aperture. 


Time(s) 


Figure  11.  Transient  behavior  of  multi-transverse- 
mode  VCSEL.  Photon  numbers  in  single  modes  and 
total  number  of  photons  in  the  cavity. 


Figure  13.  Steady-state  axial  distribution  of  intensity 
in  transverse  modes  for  VCSEL  with  6  micrometer 
oxide  aperture. 


t  =  0.44  ns 


t  =  0.48  ns 


t  =  0.60  ns 


Figure  14.  Selected  snapshots  of  transient  behavior  of  total  optical  intensity  during  the  multi- 
transverse-mode  VCSEL  switching  (cf.  Figure  11) 
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The  above  simulations  demonstrate  the  first  in  the  U.S.  commercial  software  for  high-fidelity 
multi-mode  simulation  of  lasers  (VCSELs  and  edge-emitting  lasers),  by  optical-electrieal 
coupling  of  CFD-ACE-l-  Multiphysics  Simulator  from  CFDRC  with  WIMP  Optics  from  AFRL. 

•  The  first  training  course  for  DoD  users  was  held  at  WP-AFB  on  Oct.  1-3,  2002. 

•  Two  new  tutorials  (one  for  VCSEL,  one  for  GaN  EEL)  have  been  developed  and  printed. 
This  makes  already  4  different  O’ SEMI  tutorials  available  for  users. 

•  Two  separate  licenses  of  the  O’ SEMI  Package  have  been  shipped  to  University 
California  at  Santa  Barbara  (one  to  Prof  Joachim  Piprek,  and  second  to  Prof  Larry 
Coldren).  UCSB  is  the  leading  U.S.  institution  in  the  R&D  efforts  for  new  lasers  (VCSELs) 
and  optoelectronic  materials. 

2.3  Quantum  mechanical  model  enhancement  of  electronic  and  photonic  devices 

Lead:  Cornell  University,  PI:  Edwin  Kan 

Partial-differential-equation  (PDE)  level  inclusion  of  quantum  mechanical  effects  in  nano-scale 
electronic  and  photonic  devices  is  difficult  in  direct  formulation,  mainly  due  to  nonlocality  and 
coherence  in  phase  space.  The  problem  has  been  attempted  early  in  the  development  history  of 
quantum  mechanics,  in  particular  Feynman's  effective  potential  smoothing  based  on 
Nonequilibrium  Green's  Function  (NEGF).  At  Cornell,  the  effective  potential  formalism  and 
several  other  recent  proposals  such  as  quantum  potentials  and  transmission  matrix  have  been 
compared  analytically  and  numerically  (on  PROPHET  platform).  All  formalisms  need  to  couple 
to  Poisson  equation  and  carrier  statistics  on  the  boundary  to  be  meaningful  for  device  operational 
modeling.  Self-consistent  solution  was  attempted  in  the  first  year  of  the  project,  to  understand  its 
physical  accuracy  and  numerical  properties. 

For  a  general  semiconductor  device  we  can  split  the  potential  into  the  sum  of  a  smooth  function 
and  a  discontinuous  function.  A  key  point  to  note  here  is  that  the  discontinuous  part  (or  large 
perturbation  in  terms  of  the  thermal  wavelength)  of  the  potential  is  caused  by  an  abrupt  material 
interface  and  is  known  a  priori  once  the  device  structure  is  known.  It  is  independent  of  any  self- 
consistent  calculation  required  to  calculate  the  total  potential  by  iterating  the  transport  equations 
with  the  Poisson  equation.  We  can  therefore  in  principle  calculate  the  effect  of  the  large  potential 
variations  exactly  using  the  integral  form  given  above.  In  other  words,  the  Bloch  equation  now 
becomes  Utot  =  Usm  +  Udu-  For  any  simple  1-D  barrier  given  by,  such  as  for  example  found  in  a 
MOS  transistor  or  a  quantum  well,  we  can  now  solve  for  the  Wigner  function  in  equilibrium. 
Alternatively,  we  can  try  to  Fourier  transform  the  whole  Bloch  equation  with  respect  to 
momentum  (in  effect  transforming  the  equation  to  that  for  the  density  matrix  in  the  difference 
coordinate). 

We  can  then  take  moments  on  this  equilibrium  Wigner  function  to  obtain  self-consistent 
solutions  with  the  Poisson  equation  and  inclusion  of  scattering. 
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In  the  moment  equation,  the  only  different  term  from  the  conventional  Boltzmann  transport  is  at 
the  stress  tensor: 
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This  can  be  seen  as  the  density  gradient  theory  with  a  position-dependent  correction  term,  which 
captures  the  drop  in  density  of  carriers  near  the  barrier  due  to  barrier  repulsion.  The  tunneling 
current  in  finite  barrier  and  temperature  can  thus  be  self-consistently  included.  We  are  now 
trying  to  generalize  this  formalism  to  multi-dimension  cases  in  realistic  device  structures.  The 
goal  is  to  obtain  a  tangible  PDE  form,  and  later  on,  further  reduce  to  compact  models  for 
quantum  electronic  and  photonic  devices. 

In  next  step,  we  obtained  a  full  set  of  novel  Wigner  moment  equations  shown  above.  The 
equation  is  not  a  simple  addition  of  the  quantum  correction,  but  by  moment  generation  from  an 
exact  Wigner  function  in  equilibrium  in  the  presence  of  a  large  potential  barrier  resulted  from 
material  heterojunctions.  We  have  compared  this  approach  with  the  density  gradient  theory  and 
have  found  significant  improvement  in  modeling  the  carrier  behavior  inside  the  barrier 
(tunneling  by  the  evanescent  mode).  Then,  we  formulated  various  measurable  operators, 
including  carrier  concentration  and  carrier  flux.  Details  of  the  equation  are  presented  in  the  two 
papers  on  Quantum  Hydrodynamic  Modeling  of  Semiconductor  Transport  (listed  in  References 
section),  by  Venkat  Narayanan  and  Edwin  Kan  from  Cornell  University. 

As  a  final  stage  of  corroboration  of  the  Wigner-based  moment  equation,  we  have  composed  a 
nonequilibrium  Green  function  (NEGF)  solver,  which  provides  information  in  the  dimensions  of 
space,  momentum  and  energy.  For  a  1-D  case  with  only  optical  and  acoustic  phonon  scattering, 
the  solver  can  be  put  together  in  relatively  straightforward  manner  instead  of  re-engineer 
software  from  elsewhere.  The  computational  efficiency  for  NEGF  is  however  only  good  for 
corroboration  since  any  realistic  device  structure  will  take  unrealistic  time.  The  full  Wigner 
solver  we  have  obtained  last  term  is  also  functioning  and  provide  another  benchmark.  We  have 
used  these  solvers  to  compute  the  explicit  formula  of  stress  tensors  and  energy  flow  in  the 
Wigner-based  moment  equation  and  the  density  gradient  algorithm  in  order  to  understand  the 
underlined  assumption  and  accuracy.  The  figures  below  show  the  stress  tensor  comparison  in  a 
tunneling  barrier  under  equilibrium  and  0.5  V  bias. 
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Quantum  Transport 


■  Quantum  hydrodynamic 
(QHD)  equations  are  the  same 
as  classical  ones  in  the  first 
three  moments,  and  the  only 
difference  is  at  the  closure. 

■  Add  phenomenological 
scattering  (required  at 
boundaries) 


Construct  the  density  matrix  from  NEGF, 
then  construct  the  pressure  tensor  from 
definition  (closure  for  2"^  moment) 
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Evaluate  pressure  tensor  from  the  kinetics 
(inspired  by  the  equilibrium  guess)  (QHD) 
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■  Use  NEGF  to  obtain  QHD  or 
QDD  closure  conditions 

■  Benchmark  with  existing 
approaches  such  as  density 
gradients 
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Further  simplification  to  n  only  (QDD) 
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Stress  Tensor  Comparison  (I"**  Moment  Closure) 


Figure  15.  Quantum  transport  modeling  procedure,  and  sample  results  of  the  stress  tensor  in  a 
tunneling  barrier  under  equilibrium  (0  V)  and  1  V  bias. 
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3  IMPLEMENTATION  OF  HIGH  ACCURACY  ROBUST  3-D  EM  SOLVERS 


CFDRC  worked  on  two  different  3-D  Electromagnetic  Maxwell  equation  solvers: 

•  Finite-Difference  Time-Domain  (FDTD)  Solver:  CFD-Maxwell,  and 

•  Finite-Volume  Magnetic- Vector-Potential  based  solver  (EMAG) 

They  were  modified  and  enhanced  to  be  maximally  efficient  when  used  with  Cartesian  grids,  new 
types  of  unstructured  grids,  as  well  as  with  Filament  and  Wire  models.  There  are  cases  when 
implicit  code  (EMAG)  is  better,  and  in  other  cases  explicit  FDTD  code  (CFD-Maxwell)  may  be 
more  efficient.  Examples  are: 

*  CFD-Maxwell  scales  linearly  with  the  number  of  mesh  cells,  which  is  a  big  advantage  over  the 
implicit  electromagnetic  solver  (EMAG).  For  very  large  problems  (millions  of  mesh  cells),  the 
explicit  FDTD  code  will  be  probably  the  only  available  option. 

*  CFD-Maxwell  is  intrinsically  suited  for  parallelization. 

*  EMAG  is  good  for  coupled  multi-physics  problems,  e.g.  with  thermal,  mechanical,  etc.,  as  an 
integral  part  of  the  Multi-Physics  Suite  CFD-ACE-i-.  Also,  it  should  be  more  efficient  for  stiff 
problems,  like  solving  conductive  currents  in  wires  together  with  electromagnetic  fields  in 
surrounding  insulators. 

Therefore,  we  continued  to  develop,  improve,  and  test  both  the  electromagnetic  codes  at  CFDRC, 
because  we  had  seen  a  big  potential  for  both  of  them  in  the  future  comprehensive  CAD  package  for 
Opto-Electronic  Interconnects. 

3.1  Finite-Difference  Time-Domain  (FDTD)  Solver:  CFD-Maxwell 

3.1.1  FDTD  Solver  Adaptation  and  Verification  for  RF  Interconnects 

Primary  target  was  to  obtain  sufficient  trust  in  the  reliability  and  accuracy  of  CFDRC ’s 
electromagnetic  3D  full-wave  FDTD  solver,  CFD-Maxwell.  This  was  an  essential  foundation  for 
future  developments.  Therefore,  we  performed  several  validation  studies,  described  below. 

Example  1 

The  first  validation  case  was  a  shielded  coupled  microstrip  line  problem,  as  shown  below. 


1.15  miD 


4.0  mm 

Figure  16.  The  coupled  microstrip  line.  The  strips  are  0.05  by  0.3  mm  and  separated  by  0.3  mm, 
the  dielectric  slab  is  0.25  mm  thick  and  has  a  relative  dielectric  constant  Sr  of  4.5. 
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Results  from  earlier  supercomputer  simulations  for  this  specific  case  have  been  published  by 
Becker,  Harm  and  Mittra  in:  “Time-Domain  Electromagnetic  Analysis  of  Interconnects  in  a 
Computer  Chip  Package”,  IEEE  Transactions  on  Microwave  Theory  and  Techniques,  Vol.  40, 
No.  12,  Dec.  1992.  The  results  obtained  with  CFD-Maxwell  (Figure  17  and  Figure  18)  match 
these  previously  published  numbers  with  discrepancies  on  the  order  of  1%  to  2%. 


Figure  17.  Propagation  of  an  even  mode  Gaussian  pulses  (Pulse  Width  10"”  s)  along  a  coupled 

microstrip  line. 
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Figure  18.  CFD-Maxwell  impedances  for  even  and  odd  mode  excitation  in  a  frequency  range 
between  5GHz  and  25GHz  in  comparison  with  the  results  from  Becker,  Harm  and  Mittra 

(IEEE  Trans,  on  MTT,  1992) 
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Example  2 


As  a  second  example  we  chose  a  5-line  interconnect  structure  published  in  [Beattie,  2000],  and 
shown  in  Figure  19.  It  is  a  five  line  signal  bus  over  a  ground  mesh.  For  all  signal  lines  driver 
resistance  is  50  Q,  load  capacitance  2  fF,  thickness  3.5  pm,  spacing  10  pm,  length  1000  pm  and 
width  10  pm.  Spacing  to  ground  plane  is  5  pm.  System  size  is  about  1005  pm.  Materials  are 
copper  (p=  17.5  nQm)  and  Si02  (Sr=  3.9,  pr  =1.0).  The  middle  signal  line  is  active.  In  CFD- 
Maxwell  at  that  time  we  didn’t  have  yet  available  a  lumped  capacitance  model  (it  was  developed 
later),  so  we  used  a  resistive  line  terminations  on  the  far  end. 


Figure  19.  Test  case  of  interconnects:  five-line  signal  bus  over  a  ground  mesh  (a)  and  its  transfer 
function  from  near  to  far  end  node  of  leftmost  signal  line,  calculated  with  different  equivalent 

circuit  models  in  [Beattie,  2000]. 


A  3D  computational  meshed  model  for  CFD-Maxwell,  built  automatically  using  CFD-Micromesh, 
is  shown  in  Figure  20. 
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Figure  20.  Our  meshed  model  of  the  solid  (Cu)  parts  of  the  5-line  bus  interconnects  from  Figure 
19a.  The  mesh  was  built  automatically  in  CFD-Micromesh  within  a  few  seconds. 


Figure  21.  Two  cross-sections  showing  meshed  insulator  parts  of  the  5-line  interconnect  structure 
(the  copper  parts  are  blanked  for  better  visualization).  The  full  mesh,  built  automatically  in  CFD- 

Micromesh,  has  about  46,000  cells. 
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Figure  22.  Several  selected  snapshots  from  CFD-Maxwell  transient  simulation  of  a  propagating 
Gaussian  pulse  of  Ips  width,  after  1000  time  steps  (5  ps  of  the  simulated  time),  1500  steps  (7.5  ps), 
2000  steps  (10  ps),  and  2500  steps  (12.5  ps).  The  last  snapshot  shows  a  reflected  wave  from  the  un¬ 
matched  line  termination. 

An  electric  signal  pulse  was  applied  to  the  middle  line,  from  a  voltage  source  with  internal 
resistance  50  Q,  as  illustrated  in  Figure  19a.  Two  types  of  pulses  were  applied:  Gaussian  and 
trapezoidal  (which  is  closer  to  the  shape  of  digital  signals).  Several  selected  snapshots  of  CFD- 
Maxwell  simulation  of  a  propagating  Gaussian  pulse  of  Ips  width  are  shown  in  Figure  22.  The 
plots  on  the  left  side  of  each  case  show  voltage  distribution  along  the  active  line  (upper  plot)  and 
along  the  neighbor  line  (lower  plot)  for  each  time  instant.  This  allows  for  cross-talk  analysis  in 
the  interconnect  structure.  For  example,  in  the  first  (upper-left)  snapshot  in  Figure  22,  at  time 
step  1000  (which  corresponds  to  5  ps  of  the  simulated  time)  the  cross-talk  voltage  on  the 
neighbor  line  reaches  about  20%  of  the  active  voltage. 

To  calculate  the  characteristic  impedance  of  the  line,  first  we  terminated  all  the  lines  with  a  short 
circuit  (R  =  0).  We  monitored  the  voltage  and  current  values  in  the  middle  line  at  the  near  end  of 
the  line  (v7  and  il,  at  0.1  of  the  length)  and  at  the  far  end  of  the  line  (v2  and  i2,  at  0.9  of  the 
length).  Time  dependent  waveforms  of  these  monitored  voltage  and  current  values  are  presented 
in  Figure  23.  On  the  time  scale,  1000  time  steps  corresponds  to  5  ps  of  simulated  time.  In  this 
case,  due  to  the  short  circuit  at  the  line  termination,  the  voltage  pulse  is  reflected  with  the 
opposite  sign  from  the  far  end,  and  comes  back  to  the  near  end  of  the  line. 
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Figure  23.  Voltages  (left  plot)  and  currents  (right  plot)  varying  in  time  in  the  middle  line  of  the  5- 
line  bus,  with  a  short  circuit  (R=0)  at  the  line  termination.  The  values  were  monitored  at  the  near 
end  of  the  line  (vi  and  il,  at  0.1  of  the  length  -  red  curves)  and  at  the  far  end  of  the  line  (v2  and  i2, 
at  0.9  of  the  length  -  green  curves).  1000  time  steps  corresponds  to  5  ps  of  simulated  time. 

From  the  first  simulation  we  determined  the  characteristic  impedance  of  the  line  to  be  38.7  Q. 
When  we  applied  R  =  38.7  Q  as  the  line  termination,  to  ensure  the  matching  impedance  at  the  end, 
the  signals  were  not  reflected  anymore,  as  illustrated  in  Figure  24  and  Figure  25. 


Figure  24.  Time-dependent  voltages  (left)  and  currents  (right)  in  the  middle  line  of  the  5-line  bus, 
with  a  matching  impedance  IR  =  38.7  Q )  at  the  line  termination. 


Figure  25.  Time-dependent  voltages  (left)  and  currents  (right)  in  the  middle  line  of  the  5-line  bus, 
with  a  matching  impedance  (R  -  38.7  Q )  at  the  line  termination,  with  a  trapezoidal  pulse  applied 
to  the  input  of  the  middle  line.  1000  time  steps  corresponds  to  5  ps  of  simulated  time. 
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From  the  Gaussian  pulse  simulation,  illustrated  in  Figure  24,  we  ealeulated  the  S-parameters  of  the 
middle  line  of  the  5-line  bus,  to  obtain  among  others  the  transfer  characteristic  similar  to  the  one 
presented  in  Figure  19b.  Our  calculated  S-parameters  for  the  middle  line  are  shown  in  Figure  26. 


frequency  [Hz]  (LOG) 

Figure  26.  Plot  of  S-parameters  calculated  from  the  transient  simulation  of  Gaussian  pulse. 


The  problem  with  our  calculation  of  S-parameters  as  a  function  of  frequency  from  transient  results 
(through  Fast  Fourier  Transform,  or  FFT)  is  that  the  lowest  calculated  frequency  is  proportional  to 
the  extent  of  sampling  time  of  the  time-domain  signal.  Since  our  first  transient  calculations 
featured  very  fast  signal  (Ips  wide  pulse)  and  very  short  covered  time  (total  of  20  ps),  the  lowest 
frequency  value  obtained  from  the  FFT  was  9.3E10  Hz,  or  93  GHz  (Figure  26),  which  was  high 
above  the  range  in  which  we  would  like  to  obtain  the  S-parameters  (cf  Figure  19b).  To  cover  the 
lower  frequency  range,  we  should  run  a  longer  signal  pulse,  in  a  longer  simulated  time  range. 
However,  the  time  step  in  the  FDTD  method  is  limited  by  the  smallest  mesh  cell  size  (Courant 
condition),  so  we  could  not  increase  time  step  for  this  particular  structure  (since  the  smallest  mesh 
cell  size  was  determined  by  the  smallest  physical  dimension  of  the  line,  that  is  its  thickness  of  3.5 
pm).  Therefore,  we  had  to  use  the  same  time  step,  but  just  to  run  the  simulation  for  much  more 
time  steps,  going  into  thousands.  Sample  results  of  such  simulation  of  a  longer  Gaussian  pulse  (10 
ps  wide)  are  presented  in  Figure  27a,  and  for  a  trapezoidal  pulse  10  ps  wide  -  in  Figure  27b. 
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Figure  11.  Time-dependent  voltages  (left)  and  currents  (right)  in  the  middle  line  of  the  5-line  hus, 
with  a  lOps-wide  Gaussian  pulse  (a)  and  lOps-wide  trapezoidal  pulse  (b)  applied  to  the  input  of  the 
middle  line.  On  the  time  scale,  1000  time  steps  corresponds  to  5  ps  of  simulated  time. 


In  the  cases  illustrated  in  Figure  27  we  hit  the  known  generic  problem  of  the  explicit  FDTD 
method,  called  “late  time  instability”,  which  causes  FDTD  simulations  on  not  perfectly  rectangular 
meshes  (which  is  our  case)  with  many  thousands  of  time  steps  to  diverge  and  prevents  one  from 
obtaining  a  stable  solution  for  longer  simulated  time  intervals.  Since  the  time  step  can  not  be 
increased  due  to  the  Courant  condition,  as  explained  above,  this  poses  a  serious  problem  for  longer 
times,  or  lower  frequencies,  or  EM  simulations  of  structures  containing  very  small  geometrical 
features  (which  force  a  very  small  time  step). 


An  alternative  would  be  using  CFDRC’s  implicit  electromagnetic  solver  (EMAG  module  in 
ACE-I-)  which  does  not  impose  such  strict  limits  on  the  time  step  length.  The  EMAG  solver  during 
the  NeoCAD  Program  was  enhanced  with  special  features  for  electrical  interconnect  simulations, 
like  extracting  and  displaying  time-varying  voltages  and  currents  in  metal  lines,  calculating 
characteristic  impedance,  S-parameters,  etc.  They  are  described  in  later  sections  of  this  report. 


New  Capability  in  CFD-Maxwell:  Capacitive  Line  Termination 

In  many  cases  in  integrated  circuits,  the  termination  of  an  interconnecting  line  may  have 
capacitive  character.  In  particular,  for  signal  lines  it  can  be  gate  of  a  MOSFET.  Quite  often  it  is 
the  most  appropriate  to  model  the  interconnect  transmission  lines  as  driven  by  a  resistive  voltage 
source  at  the  input,  and  terminated  by  load  capacitances,  as  illustrated  in  Figure  28. 
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C,  termination 


R,  source 


Figure  28.  Example  of  five  parallel  signal  lines  driven  by  resistive  voltage  sources  and  terminated 

by  load  capacitances. 


The  characteristic  impedance  of  the  transmission  line  determines  the  current  that  can  be  supplied 
to  charge  capacitive  load  Cl.  From  the  load’s  point  of  the  view,  the  line  behaves  as  a  resistance 
with  value  Zq.  The  transient  response  at  the  capacitor  node,  therefore,  displays  a  time  constant 
ZqCl.  This  is  illustrated  in  Figure  29,  which  shows  the  simulated  analytical  transient  response  of 
a  transmission  line  with  a  characteristic  impedance  of  50  Q  loaded  by  a  capacitance  of  2  pF.  The 
response  shows  how  the  output  rises  to  its  final  value  with  a  time-constant  of  100  psec  (=50 
Q  X  2  pF)  after  a  delay  equal  to  the  time-of-flight  of  the  line. 


Figure  29.  Theoretical  transient  response  of  a  transmission  line  with  a  characteristic  impedance  of 
50  Q,  loaded  by  a  capacitance  of  2  pF,  to  an  applied  step  voltage  of  5V. 

To  verify  the  newly  implemented  capacitive-termination  capability  in  CFD-Maxwell,  we  have 
performed  simulations  of  an  applied  5V  step  voltage  to  calculate  transient  response  of  a 
transmission  line  loaded  by  a  capacitance  and  by  a  matching  impedance  (resistance).  The  obtained 
waveforms  of  input  voltage  (Vin),  output  voltage  (Vout),  input  current  (lin),  and  output  current 
(lout)  for  the  line  with  capacitive  termination  are  shown  in  Figure  30.  Our  FDTD  results  of  voltage 
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waveforms  are  very  similar  to  the  expected  theoretical  results  from  Figure  29.  This  confirms 
proper  implementation  of  the  new  capacitive  line  termination  model  in  CFD-Maxwell. 


Figure  30.  CFD-Maxwell  results  for  a  transmission  line  with  capacitive  termination; 
time  dependent  voltages  (left)  and  currents  (right). 


3.1.2  FDTD  Solver  Applications  for  Optical  Modeling 

The  tests  reported  above  focused  on  examples  verifying  the  validity  of  solutions  for  basic,  high 
frequency  electronic  building  blocks  (transmission  lines.)  Equally  important,  however  more 
complicated,  is  the  validation  of  pure  wave  propagation  problems  since  in  those  cases  not  single 
quantities  (such  as  impedance)  or  one-dimensional  functions  (such  impedance  as  a  function  of 
frequency)  are  the  results  of  simulations,  but  in  many  cases  2D  or  3D  wave  propagation  patterns 
and/or  energy  propagation  distributions.  Such  solutions  are  inherently  sensitive  to  mesh 
imperfections  and  numerical  stability  issues  have  to  be  observed  carefully. 

Moreover,  beside  obvious  cases  of  instability,  FDTD  simulations  can  yield  reasonably  appearing 
results  that  are  inaccurate.  Rules  have  been  established  for  acceptable  mesh  configurations  and 
aspect  ratio  changes,  including  guidelines  concerning  the  mixed  application  of  structured  and 
unstructured  meshes,  since  such  rules  are  not  yet  incorporated  in  intelligent  agents  that  would 
control  the  meshing  process  (Figure  31  through  Figure  33).  Simulation  results  must  be  trusted 
sufficiently  that  they  can  form  the  basis  for  devices  to  be  fabricated. 

The  application  of  shape-conformal  meshes  is  imperative  for  precise  wave  scattering  solutions. 
This  capability  is  one  of  the  key  features  provided  by  CFD-Maxwell  and  was  utilized  for  the  OE- 
IDEAS  developments.  In  order  to  confirm  proper  operation  of  the  code,  examples  with  analytical 
solutions  are  invaluable.  EM  wave  scattering  on  dielectric  and  perfectly  electrically  conducting 
cylinders  are  two  of  the  few  known  cases  for  which  analytical  solution  involving  infinite  series  of 
Bessel  and  Hankel  functions  and  satisfying  Maxwell’s  equations  can  be  found. 
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Figure  31.  TM  wave  propagation  test  in  2D  domain,  A,  /Ax  =12  [Ax  =cell  edge  length],  solver  in  total 
field  formulation  mode,  E^(x,y)  plotted  at  arbitrary  time.  Left:  unstructured  mesh,  left  and  right 
boundary  ABC,  strong  distortion  of  wave  fronts;  Right:  unstructured  mesh  +  Cartesian  buffer  zone, 
left  and  right  boundary:  extrapolation,  wave  fronts  are  preserved. 


Figure  32.  TM  wave  propagation  test  in  narrowing  2D  Figure  33.  Initial  TM  wave  scattering 
domain.  A,  /Ax  »  20,  wave  fronts  are  preserved.  test  on  a  dielectric  square  with  Sr=4, 

solver  in  total  field  formulation  mode. 


Cylinders  provide,  therefore,  an  important  initial  object  to  study  the  accuracy  of  scattering 
solutions.  Ultimately,  the  systems  of  interest  may  contain  a  large  number  of  scatterers.  It  is, 
therefore,  of  value  to  find  the  minimal  required  mesh  density  for  a  desired  level  of  accuracy. 

Using  a  computer  algebra  system,  it  is  possible  to  map  time-dependent  analytical  wave  solution 
for  complete  2D  spaces,  including  inside  materials,  and  compare  these  results  to  FDTD  data 
(Figure  34  through  Figure  36).  Due  to  the  use  of  shape-conformal  meshes,  accuracies  on  the 
order  of  2%  for  the  steady  state  field  values  have  been  found  for  a  number  of  cases  involving 
different  wavelength-to-scatterer  size  ratios  and  permittivities.  Gradually,  using  larger  and  more 
complex  examples,  strong  confidence  in  the  solver  has  been  established.  A  true  3D  scattering 
problem  with  over  10*’  cells  is  shown  in  Figure  37. 
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Figure  34.  Scattering  of  a  plane  electromagnetic  wave  on  an  infinite  dielectric  cylinder  with  £>=4, 
diameter  d  (normal  incident,  TM^  polarization,  +x-propagation);  /L4?=1.5;  (solver  in  scattered  field 
mode,  scattered  field  plotted)  left:  FDTD  simulation,  contour  plot  of  (r,^2n  t/T)\  at  constant, 
random  phase  angle  rafter  steady  state  is  reached;  middle:  2D  mapping  of  the  analytical  solution 
for  \^(r,^2n  t/T)\  at  same  phase  angle;  right:  surface  plot  of  analytical  solution  for  ^(r,^2n  t:/T)\ 


Figure  35.  Left:  Comparison  of  analytical  and  FDTD  solution  for  F’/(x)|y=o;  middle:  FDTD 
simulation,  contour  plot  of  total  field  \^(r,^,2n  t/T)\  ;  right:  2D  mapping  of  the  analytical  solution 

for  total  field  \^(r,^2n  t/T)\ 


Figure  36.  Scattering  of  a  plane  electromagnetic  wave  on  an  infinite  PEC  cylinder  with  diameter  d 
(normal  incident,  TE*  polarization,  +x-propagation);  A/d=\.S‘,  (solver  in  scattered  field  mode, 
scattered  field  plotted)  left:  FDTD  simulation,  contour  plot  of  \sB(r,^f27t  t/T)\  at  constant,  random 
phase  angle  rafter  steady  state  is  reached;  middle:  comparison  of  analytical  and  FDTD  solution  for 
sBj(x)\y^  at  same  phase  angle;  right:  2D  mapping  of  the  analytical  solution  for  \sB(r,^27t  t/T)\ 


27 


Figure  37.  3D  FDTD  simulation  of  the  scattering  of  a  TM*  plane  wave  packet  (propagating  in  +z 
direction  coming  from  the  bottom  of  the  image,  gray  area)  through  a  2x2  square  aperture  array 
(each  aperture’s  size  is  Ax/L);  walls:  contour  plot  of  (x,y,z,t)\,  volume:  iso-surfaces  of 
(x,y,z,t)\=QA\ Im  illustrating  the  beam-lets  exiting  the  apertures,  notice  the  formation  of  the  main 

diffraction  maxima  in  the  middle. 


3.2  Finite-Volume  Magnetic- Vector-Potential  based  3D  Solver;  EMAG 
3.2.1  EMAG  Solver  Modifications  (Voltage,  Current,  Impedance,  S-Parameters) 

Modifications  to  the  EMAG  module  in  CFD-ACE+  were  aeeomplished,  allowing  for  the 
monitoring  of  voltage  and  current  values  in  three-dimensional  simulations.  The  voltage  and 
current  ealeulations,  when  aetivated,  are  output  at  the  end  of  eaeh  time  step  in  a  transient 
simulation  and  at  the  end  of  a  steady  state  simulation.  The  output  voltage  and  eurrent  values 
may  then  be  used  for  impedance  and  S-parameter  ealeulations. 

The  voltage  and  current  output  may  be  activated  through  DTE  file  updates,  aetivating  the  voltage 
and  current  output  capabilities.  Voltage  is  calculated  by  integrating  the  total  eleetrie  field  along 
a  line  eonnecting  the  users  desired  measurement  point  and  a  reference  ground  point.  The  current 
ean  be  calculated  in  two  ways.  The  current  loop  technique  integrates  the  magnetic  field  vector 
along  a  loop  around  the  region  of  interest.  The  current  face  technique  integrates  the  total  eurrent 
density  vector  across  an  interface.  The  impedance  and  S-parameter  ealeulations  are  ealculated 
using  CFD-Maxwell’s  post  processing  executables.  The  teehniques  are  sehematieally  illustrated 
(along  with  their  basie  equations)  in  Figure  38. 
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Figure  38.  Schematic  illustrating  voltage  lines,  current  loops,  and  current  faces  (basic  equations  for 

voltage,  current,  and  impedance  are  included). 

3.2.2  Coupling  of  3D  EMAG  Simulator  with  SPICE 

Under  the  NeoCAD  funding,  the  new  numerical  procedures  allowing  for  direct  coupling  of 
CFDRC’s  3D  Electromagnetic  Solver  (EMAG)  with  the  SPICE  circuit  simulator  were 
developed,  implemented,  and  tested.  The  coupling  mechanism  is  illustrated  below. 
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(integrated  from  BC  faces) 


Figure  39.  Illustration  of  the  coupling  mechanism  of  CFDRC  3D  EMAG  simulator  with  SPICE. 

Coupling  Algorithm 


The  coupling  of  SPICE3  (from  UC  Berkeley)  with  CFD-ACE+  Multiphysics  Simulator  was 
achieved  through  a  two-level  Newton  algorithm.  All  ACE  models  become  a  new  type  of 
element  in  the  SPICE  netlist.  Linearized  models  of  ACE  devices  are  computed  by  the 
CFD-AceSpice  program  (described  below)  at  each  step.  The  Jacobian  is  obtained  by  numerical 
perturbations  of  the  ACE  solution,  or  from  previous  steps. 
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Figure  40.  Flow-chart  of  the  SPICE  to  ACE+EMAG  coupling  algorithm. 


CFD-AceSpice  program  or  “controller” 

CFD-AceSpice  is  the  name  given  to  the  interface  between  the  CFD-ACE+  and  SPICES 
programs.  It  starts  instances  of  both  solvers,  co-ordinates  information  exchange  between  the  two 
solvers,  and  directs  the  progress  of  the  simulation  based  on  the  convergence  criteria,  tolerances 
etc.  specified  in  the  input  circuit  (.cir)  file. 

The  simulation  proceeds  as  follows: 
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•  CFD-AceSpice  starts  SPICES  and  ACE  solvers.  Only  one  instance  of  SPICES  and  one  ACE 
instance  per  A-type  device  is  used  in  the  SPICE  netlist. 

•  The  SPICES  process  is  always  started  on  the  local  machine,  while  ACE  programs  can  be 
started  locally  or  on  a  remote  server. 

•  For  every  time  step  in  a  transient  analysis  (.TRAN),  or  DC  operating  point  in  a  DC  sweep 
point  analysis  (.DC),  a  Newton-Raphson  iteration  loop  is  run  until  convergence. 

This  mechanism  allows  to  run  simulations  of  complex  circuits  described  by  a  SPICE  netlist 
directly  coupled  with  a  SD  full-wave  EM  solution  of  critical  parts  (e.g.  interconnects,  vias, 
connectors,  etc.)  of  an  RF  integrated  circuits  or  device.  This  idea  is  schematically  illustrated 
below,  where  driver  and  termination  are  represented  by  SPICE  circuit,  while  the  transmission 
line  interconnects  are  modeled  by  a  3D  EMAG  model. 


Figure  41.  An  example  of  a  3D  EMAG  model  of  a  transmission  line  with  SPICE  modeling  the 

source  (driver)  and  the  line  termination. 


This  new  EMAG-SPICE  coupling  software  was  offered  to  partners  in  the  NeoCAD  program  for 
testing  and  evaluation. 
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3.2.3  3D  EMAG  Simulation  of  a  Parasitic  Hanging  Stub  at  Vias  in  RF  PCB  Interconnects 

The  phenomenon:  the  transfer  function  of  an  interconnect  trace  on  a  printed-circuit  board  (PCB) 
backplane  of  an  electrical  system  can  be  severely  dumped  at  some  frequency  due  to  a  floating  via 
stub  which  is  part  of  a  through-hole  via.  Most  fabricated  PCB  transition  types  leave  an  open  via 
stub  (see  Figure  42).  The  problems  caused  by  this  parasitic  stubs  were  shown  on  experimental 
data  of  PCB  transfer  function,  presented  at  2003  IEEE  International  Solid-State  Circuits  (ISSC) 
Conference  in  San  Francisco,  Feb  2003  (Figure  43). 

This  problem  was  analyzed  in  detail  during  the  OE-IDEAS  project  by  the  University  of 
Delaware  group.  A  detailed  description  of  different  approaches  in  approximate  modeling  of  the 
via-stub  structure  using  HSpice  models  was  presented  in  one  of  the  quarterly  reports. 


Trace  on  layer  #1 


n 


VIA 

Trace  on  layer  #2 


Via  through  hole 


< —  Stub  near  end 
STUB 

< —  Stub  far  end 


Figure  42.  Schematic  side-view  drawing  of  a  PCB  via  with  hanging  stub. 


Maximum  Stub  vs  Minimum  Stub  (Standard  0.D18"  PTH,  12in  FR-4,  0.240"  Thick) 


Figure  43.  Measured  transfer  characteristics  (S21)  of  via  with  hanging  stub  [ISSC  2003]. 
The  resonance  dip  in  the  blue  curve  is  caused  by  a  long  (240  mils)  hanging  via. 
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In  this  section,  we  summarize  a  full-wave  3D  electromagnetic  modeling  of  the  PCB  via  with  and 
without  parasitic- stub  structure  using  the  EMAG  tool  at  CFDRC. 

First,  a  layout  of  the  PCB  structure  with  two  transmission  lines  and  via  (similar  to  Figure  42)  was 
drawn  in  CFD-Micromesh  layout  editor;  see  Figure  44. 


Figure  44.  Layout  of  two  transmission  lines  and  via  (similar  to  Fig.  1),  drawn  in  Micromesh. 

From  this  layout.  Micromesh  built  automatically  a  3D  computational  model  and  3D  mesh  from 
the  technology  data  of  the  PCB,  signal  traces,  and  via,  provided  by  University  of  Delaware.  A 
close-up  view  of  a  fragment  of  the  meshed  model  is  shown  in  Figure  45. 
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Figure  45.  A  fragment  of  the  3D  meshed  computational  model  generated  by  Micromesh. 

There  is  a  ground  plane  between  signal  traces  1  and  2,  not  shown  in  Figure  42. 

A  Gaussian  pulse  of  100  ps  duration  was  applied  to  the  beginning  of  the  top  transmission  line, 
and  transient  simulation  of  signal  propagation  was  performed  using  the  CFDRC  EMAG  solver. 
The  simulation  took  about  1  hour  on  a  2GHz  Pentium  PC.  The  resulting  voltage  and  current 
waveforms  were  used  to  calculate  the  characteristic  impedance  and  S-parameters  of  the  total 
transmission  line,  including  the  via. 

Next,  we  built  another  model  of  the  structure,  this  time  with  the  hanging  stub  below  the  via.  A 
view  of  this  model  with  sample  results  of  the  electric  potential  distribution  at  t  =  1  ns,  are  shown 
in  the  figure  below.  This  3D  model  was  much  larger,  and  the  simulation  with  the  hanging  stub  of 
240  mils  took  about  2.5  hours  on  a  2GHz  Pentium  PC. 
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Figure  46.  The  model  of  via  with  the  hanging  stub  (240  mils),  with  overlaid  results  of  the  electric 
potential  distribution  at  t  =  1  ns,  calculated  by  transient  EMAG  simulation.  The  ground  plane  and 
dielectric  materials  are  not  shown  here  for  better  visibility  of  the  transmission  lines. 

Also  for  the  second  structure,  with  the  hanging  parasitic  stub,  the  resulting  voltage  and  current 
waveforms  were  used  to  calculate  S-parameters  of  the  entire  transmission  line.  Comparison  of 
the  computed  S2 1  transfer  characteristics  for  the  cases  without-stub  and  with-stub  are  shown  in 
Figure  47.  The  character  of  the  simulated  transfer  characteristics  is  similar  to  the  measured  ones 
(Figure  43),  however  we  still  need  to  investigate  why  our  EMAG  simulator  calculates  S21  values 
greater  than  1,  which  are  not  physical.  The  correctness  of  detailed  physical  results  of  EMAG 
solver  is  being  analyzed  now  at  CFDRC. 
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Figure  47.  Comparison  of  the  S21  transfer  characteristics  for  the  cases  without-stub  and  with-stub, 
calculated  form  the  full-wave  transient  results  of  EMAG  simulator  from  CFDRC. 


4  IMPLEMENTATION  OF  FILAMENT  AND  WIRE  MODELS 

In  this  task,  the  proposed  filament  and  wire  models  were  being  implemented  into  EMAG,  the  3D 
electromagnetic  module  of  CFD-ACE+.  The  electric  potential  is  solved  in  each  wire/filament,  and 
their  contribution  to  the  magnetic  field  is  fed  into  the  magnetic  solver  as  source. 

4.1  Local  Subcell  Model  of  Filament  in  Electromagnetic  Module 

One  of  the  challenges  facing  the  grid-based  numerical  model  is  the  distance  scale  over  which 
physical  processes  or  material  properties  must  be  resolved  ranging  over  several  orders  of 
magnitude.  One  has  to  either 

(a)  use  a  mesh  with  globally  variable  array  of  cells  to  precisely  model  the  shape  and  feature 
of  a  structure;  or 

(b)  use  a  more  uniform  mesh  and  approximate  the  physical  properties  of  the  fine  spatial 
details  by  somehow  binding  them  into  the  local  mesh  cells. 

The  filament  model  belongs  to  the  second  type.  The  following  is  to  provide  insight  on  how  to 
model  the  contribution  of  fine  filament  to  the  much  larger  grid  system. 

Contribution  of  Conducting  Wire  to  the  Neighbor  Cells 

To  understand  the  interaction  of  a  current  carrying  wire  or  a  conductive  wire  (not  active  current) 
with  the  surrounding  meshing,  we  can  start  with  Maxwell  equation 
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VxH  =  s  — 
dt 


+  J 


where  J  is  the  total  current  density,  and  is  equal  to  oE  where  E  is 


E  =  -V(l)- 


dA 

dt 


where  A  as  magnetic  vector  potential. 

For  a  current  carrying  wire,  the  contribution  is  a  source  to  the  local  cell  with:  J ^  =  -(N<j) . 


For  a  wire  not  carrying  current,  it  could  behave  as  an  inductor  and  an  capacitor  represented  by 
dA  A 

<7  —  and  s  — T- .  There  is  no  obvious  solution  as  how  the  A  field  in  the  wire  and  the  A  field  in 
dt  g,i 

the  surrounding  medium  interacts. 

Role  as  an  Inductor 

The  contribution  of  inductive  property  of  a  conductive  wire  to  the  surrounding  medium  is 
dA 

represented  by  cr  — .  It  is  dependent  on 
dt 

(a)  conductivity  of  the  wire,  which  is  typically  very  high,  6.6  x  lO^Qm  for  copper;  and 

(b)  frequency  range  — . 

dt 


Apparently  for  high  frequencies  the  inductive  contribution  can  be  significant. 
Role  as  a  Capacitor 


The  contribution  is  e  — —  and  has  more  to  do  with  the  traveling  speed  of  light.  It  is  significant 

dE 

only  if  the  wave  length  is  compatible  with  the  wire  length  and  when  one  is  interested  in  the 
radiation  of  RF  for  the  connectors.  If  the  interconnect  is  1  mm  the  frequency  has  to  be  300  GHz 
before  this  term  becomes  significant.  From  this  point  of  view,  we  can  negate  the  contribution  of 
capacitor  of  the  wire. 

How  to  Solve  Vector  Magnetic  Potential  in  a  Filament 

Solving  magnetic  field  inside  the  conductive  wire  is  not  trivial.  All  boundaries  have  to  exchange 
information  with  it  neighboring  cells.  This  exchange  is  diffusive  (not  convective  as  in  flow  or 
heat).  In  the  convective  exchange  a  source  can  be  applied  as 
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a 


^wire 


For  diffusive  links  it  is  very  difficult  to  set  an  accurate  link.  A  good  approach  is  to  gain  some 
insight  by  studying  detailed  behavior  of  a  wire  in  a  constant  time  varying  magnetic  field. 

Computational  Model 


The  computational  model  and  grid  are  shown  in  Figure  48.  The  wire  has  a  50  pm  radius  and  a 
conductivity  of  ct  =  6.66  x  10‘’/Qm.  The  side  wall  is  symmetric  and  top  wall  is  applied  a  time 
varying  vector  potential  of 

=  sincit 


and  the  bottom  wall 


=  0.5sin(yt 


(a)  model 


(b)  grid 


Figure  48.  Computational  Model  of  a  Conductive  Wire  in  Time  Varying  Magnetic  Field 


Figure  49  shows  the  distribution  of  vector  potential  due  to  the  existence  of  the  wire.  At  low 
frequency  the  wire  has  no  affect  on  the  neighboring  media.  This  suggests  that  we  do  not  need  to 

solve  for  A  in  the  filament  and  whole  cells  can  be  treated  as  homogeneous.  The  A  value  in  the 
filament  is  the  local  value. 
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With  increasing  frequency,  the  field  around  the  wire  starts  to  change.  What  is  significant  is  that 

at  IMHz  the  A  field  in  the  wire  is  uniform  and  has  value  of  zero.  Filaments  with  high 
conductivity  and  in  high  frequency  conditions  behave  as  “perfect”  conductors.  The  same  case 
was  run  with  the  wire  removed  from  the  domain,  and  the  boundary  specified  as  perfect 
conductor.  The  result  shows  favorable  agreement  with  those  where  the  wire  is  included  (see 
Figure  50). 


Az  r 


(a)  f  =  1  kHz 
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(b)  f  =  10  kHz 


(c)  f  =  100  kHz  (d)  f-  1  MHz 

Figure  49.  Distribution  of  a  Magnetic  Vector  Potential  Around  a  Conducting  Wire 
under  a  Time  Varying  Magnetic  Field 
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(a)  perfect  conductor  lossy  conductor 

Figure  50.  Comparison  of  Field  with  Wire  Treated  as  (a)  Perfect  Conductor  and 

(b)  Lossy  Conductor 


Subcell  Modeling  of  Filament 

For  the  above  analysis  we  can  conclude  that  there  is  no  need  to  solve  A  field  inside  the  fdament. 
At  lower  frequencies  it  will  be  the  same  as  the  neighboring  media;  at  higher  frequencies  it  will 

be  a  perfect  conductor  with  A  =  0.0. 

Questions  remained  however  as  to  how  we  treat  the  “source”  of  the  filament  to  the  neighboring 
medium.  We  proposed  to  use  the  linear  properties  of  the  Maxwell  equation  and  assume  that  the 
total  field  is  the  summation 


total  ^media  ^filament 


+  A , 


Our  task  was  to  express  A  filament  analytical  solution  with  consistent  boundary  condition. 

While  Ayy^g^fa  is  ^  solution  without  existence  of  the  filament. 


If  we  define  a  local  coordinate  with  x'y'  in  the  cross-section  of  the  filament  and  z'  is  along  the 
filament  length,  the  analytical  solution  is 

V^A^  =0.0 

for  the  cylindrical  coordinate 


\  d  (  dA^\ 
r  dr  \  dr  J 


=  0.0 


A^  =  -A^gii  at  r  =  n,  a  =  radius  of  the  filament 
A^  =0  at  r=b 
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The  solution  to  the  above  equation  is 


=  Cl  +  C2  In  r 


by  substituting  the  boundary  conditions  we  have 


A=- 


Inb 


In  r 


Inb- In  r 
^  Inb  -  In 
y^lnb- In  a  y 


cell 


In  a -  Inb 


^cell 


^cell  =  - 


In^  jr 


A 


cell 


a 

A^gll  is  the  local  cell  solution  for  CFD-ACE+.  The  above  A^  represents  the  contribution  of  a 
conducting  wire  under  high  frequency. 


The  above  solution  will  ensure  that  for  r  =  a,  =  0.0.  We  can  express  the  total  solution  in 

the  following  graphic. 


total  solution 


solution  without  wire 


Az  r 


wire  solution  with 
given  B.C. 
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Filament  Wire  as  a  Source  for  Electromagnetic  Module 

In  the  following,  we  will  show  some  results  from  the  implementation  of  the  derivation:  fdament 
wire  behaves  as  a  source  to  electromagnetic  module. 

Here  the  filament  will  provide  a  source  in  the  form  of 

Jg  =(tE  =  -cN  (/> 

Jg  can  be  obtained  by  either  solving  for  electrical  DC/ AC  field  in  the  filament,  or  by  prescribed 
J g  vector. 

In  DC/AC  field  approach,  electrical  module  has  to  be  activated  and  we  are  solving: 


42 


V- J  =  0.0^VoV^^  =  0.0 


when  (/)  is  obtained,  the  electrical  field  E  =  -'V(^  can  be  calculated  and  hence  J  =  aE  . 

On  the  other  hand,  can  be  specified  as: 

j  _^j  I  ^axis  ^coil 
l^axwl  Ycoil\ 

where  is  the  vector  along  the  coil  axis  and  r^gn  is  any  point  on  the  axis  to  the  coil. 

Straight  Wire 

The  model  is  shown  in  Figure  51.  The  wire  length  is  1  and  the  voltage  difference  is  1,  so  that 
E  =  \z  .  The  filament  has  a  radius  of  0.1  with  conductivity  of  10^/Qm.  The  analytical  solution  is: 

total  I  =  A  -  J  =  Tzr'^  ■  <7\e\  =  Tir'^  x  1 0^ 

and  the  B  field  is  then 

jul  4;r  X  10”^  *;r(0.l)^  X  10^  2;tx10~^ 

B  = - = - = - 

Ittt  Itit  r 

We  considered  only  a  quarter  of  the  computational  domain.  Three  grids  were  considered: 
50x50x10,  25x25x10,  and  15x15x10.  The  grid  sizes  were  0.1,  0.2,  0.33,  respectively. 
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>  Z 


Figure  51.  Computational  Model  for  a  Straight  Wire 


Figure  52  shows  the  potential  distribution  along  the  wire  and  the  magnetic  vector  around  the 
wire. 


inii 


n 


(b)  25x25x10  grid  (c)  50x50x10  grid 

Figure  52.  Magnetic  Vector  Around  a  Current  Carrying  Wire 
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Comparison  with  analytical  solutions  are  shown  in  Figure  53.  Note,  exeept  for  the  first  point 
which  contains  the  filament,  the  solutions  mateh  with  the  analytical  solutions  exceptionally  well. 
The  inaeeuraey  in  the  first  point  is  due  to: 

1 .  It  is  the  average  solution  over  the  first  eell.  This  eell  eontains  the  wire  contribution, 
but  it  could  not  resolve  the  detailed  wire  eontribution. 

2.  Magnetie  vector  is  V  x  A,  which  depends  on  the  local  solution  A.  If  A  is  off,  so  is 
VxA. 

In  conclusion,  the  filament  model  can  provide  an  exeellent  alternative  for  eomplex  geometry  of 
wire  without  sacrifieing  aceuraey  away  from  the  wire. 


Figure  53.  Comparison  of  Magnetic  Field 


4.2  Modeling  of  Interaction  between  Filament  and  3D  EM  Field 

A  3D  model  was  built  and  run  to  validate  the  magnetie  source  through  eleetrie  eurrent  in  a  straight 
wire.  The  purpose  was  to  give  a  general  idea  on  how  the  two  are  coupled.  The  sample 
configuration  is  shown  in  Figure  54  below,  along  with  the  ealeulated  resulting  magnetic  field 
vectors. 
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Figure  54.  a)  Distribution  of  voltage  due  to  conduction  current  in  a  wire,  b)  Magnetic  field  vectors 

resulting  from  the  conduction  current  in  a  wire. 


The  idea  was  to  use  relatively  large  grid,  yet  at  the  same  time  to  be  able  to  resolve  the  sub-cell 
scale  physics  due  to  one  or  more  conducting  wires. 

Inductive  Current  in  a  Wire  due  to  Time- Varying  Magnetic  Field 

A  more  general  case  was  tested  with  fine  scale  model  for  detailed  inductive  current  of  three 
conductors  near  one  conductive  wire,  see  the  figure  below.  This  kind  of  structure  is  quite  typical 
for  multi-layered  VLSI  interconnects. 


Figure  55.  A  geometrical  setup  to  investigate  inductive  influence  of  one  wire  on  its  neighbors. 

AC  current  was  applied  to  the  top  wire  (circular  cross-section),  while  we  wanted  to  investigate  the 
induced  current  in  its  neighboring  conductors  (square  cross-section).  We  hoped  it  would  shed 
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light  on  how  to  model  the  inductive  effect  of  a  fine  wire,  or  multiple  wires,  in  a  large  control 
volume.  The  calculated  current  vectors  induced  by  the  current  flow  in  the  top  wire  are  shown  in  the 
figure  below. 


I 


Figure  56.  Computed  current  vectors  induced  by  the  current  flow  in  the  top  wire. 


EMAG  Examples 

The  multi-filament  model  can  be  used  for  all  kind  of  diffusive  physics,  including  the 
electromagnetic  vector  potential  solved  by  the  EMAG  module  of  CFD-ACE-i-.  The  new  multi¬ 
filament  diffusion  model  was  tested  against  magnetic  field  problems.  Electromagnetic  field  in  a 
square  background  with  3  thin  wires  (magnetic  field  is  generated  by  AC  current  in  one  of  the  thin 
wires)  was  simulated  with: 

1)  A  multi-domain  direct  grid  approach; 

2)  Single-filament  model  with  fine  grid; 

3)  Multi-filament  model  with  fine  grid; 

4)  Single- filament  model  with  coarse  grid; 

5)  Multi-filament  model  with  coarse  grid; 

The  grids  used  in  the  computation  were  the  same  as  in  the  test  heat  transfer  problems.  With  fine 
grid,  each  cell  has  at  most  one  filament.  With  coarse  grid,  all  3  filaments  are  in  the  same  cell.  The 
AC  fi'equency  is  1000  Hz,  the  wire  cross-section  is  0.04  x  0.04  m^,  and  the  original  current  density 
is  lOOOOOOA/m^. 

Figure  57  shows  the  typical  magnetic  field  calculated  using  the  direct  grid  approach.  Figure  58 
compares  the  magnetic  fields  along  x-axis  for  all  the  5  cases.  Table  1  compares  the  induced  current 
for  all  the  cases. 
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Figure  57.  Typical  magnetic  field  induced  by  one  of  3  wires  in  a  transmission  line. 


x(m) 


Figure  58.  Comparison  of  simulation  results  of  magnetic  fields  along  x-axis  for  5  cases  of  single- 
and  multi-filament  models,  with  fine  and  coarse  grids. 


48 


Table  1.  Induced  current  (A/m^)  for  5  cases  of  single-  and  multi-filament  models, 

with  fine  and  coarse  grids 


wire 

(source) 

2"^  wire 

3^^  wire 

Grid  solution 

485672 

205044 

146204 

Single-filament  model  with  fine  grid 

508451 

194500 

136805 

Multi-filament  model  with  fine  grid 

465474 

203200 

145794 

Single-filament  model  with  coarse  grid 

573793 

232000 

131679 

Multi-filament  model  with  coarse  grid 

425276 

210498 

210498 

Simulation  results  show  that  the  multi-filament  solution  is  close  to  the  grid  solution,  and  the 
single-filament  model  solution  when  the  grid  is  fine  enough.  For  very  coarse  grid,  when  multi- 
filaments  occupy  the  same  background  grid  cell,  the  multi-filament  model  solution  has  the 
similar  overall  accuracy  as  single-filament  model  in  the  far  field,  and  both  are  a  little  bit  off  the 
grid  solution.  However,  the  multi-filament  solution  can  reproduce  the  local  variations  of  the 
filaments  within  the  same  cell,  which  cannot  be  resolved  with  single-filament  model. 


5  AUTOMATIC  3D  MODELING  AND  MESHING  FROM  IC  LAYOUTS: 

CFD-MICROMESH 

During  the  DARPA/MTO  TEAM/NeoCAD  Summer  Project  Review,  Aug  27-29,  2002,  in 
Providence,  RI,  several  teams  indicated  the  need  for,  or  planned  the  development  of,  grid 
generation  tools.  The  round  table  discussion,  organized  by  DARPA  PMs  and  DoD  PMs,  with 
NeoCAD  Pis,  indicated  that  an  automated  geometry  modeling  and  meshing  tool,  working 
directly  from  electronic  layouts,  could  be  of  great  benefit  to  most  of  the  teams.  A  common  tool 
could  save  development  costs,  reduce  development  time,  enable  communication  between 
different  tools,  and  could  simplify  the  design  task  of  other  DARPA  programs,  e.g.  TEAM, 
MEMS,  WBGS,  and  others.  In  the  previous  couple  of  years,  DARPA  Composite-CAD  program 
partially  sponsored  CFDRC  for  the  development  of  the  Micromesh  code  for  MEMS  modeling.  It 
was  expected  that  with  a  possible  modest  enhancement  effort  the  Micromesh  could  be  used  by 
all  NeoCAD  teams  as  well  as  DoD  agencies.  A  document  summarizing  those  discussions  has 
been  prepared  by  John  W.  Rockway,  SSC-SD,  and  Andrzej  Przekwas,  CFDRC,  and  distributed 
to  all  the  interested  NeoCAD  participants.  After  that,  we  discussed  internally  at  CFDRC  the 
challenge  of  providing  CFD-ACE+  and  Micromesh  to  the  DARPA  teams.  After  the  NeoCAD 
Review  in  Providence,  and  subsequent  exchange  of  emails  within  the  interested  group  of 
NeoCAD  members,  CFDRC  made  substantial  progress  in  resolving  both  the  issues. 

As  far  as  Micromesh  is  concerned,  CFDRC  made  the  following  progress: 

•  improved  the  automating  meshing  quality  from  imported  GDSII  and  CIF  layouts, 

•  implemented  better  handling  of  large  layouts. 
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•  improved  and  simplified  API  to  DTP  mesh  fdes  to  extract  FEM  (volume)  and  BEM 
(surface)  type  meshes, 

•  a  new,  detailed  tutorial  how  to  use  Micromesh  effectively  for  NeoCAD  benchmark  example 
cases,  was  prepared  and  distributed  to  all  interested  NeoCAD  members  (see  below  selected 
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Figure  59.  Example  from  the  new  Micromesh  Tutorial,  the  Spiral  Transformer  (built  of  two  spiral 
inductors):  two  separate  GDSII  layouts  imported  into  Micromesh  from  coupled.gds  file  provided  by 

NCSU  as  a  NeoCAD  benchmark. 
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Figure  60.  3D  solid  model  of  the  NCSU  Spiral  Transformer,  built  automatically  from  layout  in 

Micromesh. 


Figure  61.  3D  prism-hex  mesh  built  automatically  from  layout  in  Micromesh. 
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We  obtained  excellent  feedback  from  Andreas  Cangellaris’  University  of  Illinois  group,  and 
several  other  NeoCAD  Pis,  in  particular  the  University  of  Washington  team,  i.e.,  University  of 
Washington  (Vikram  Jandhyala)  and  Lawrence  Livermore  National  Laboratory  (Dan  White). 

In  January  2003,  CFDRC  received  from  John  Rockway  an  interesting  report  “Comparison  of 
Tools  for  DARPA  CoSMoS  mesh  generation”  containing  collection  of  observations  of 
meshing  tools  applied  to  NeoCAD  problems.  The  observations  were  from  University  of 
Washington  and  Lawrence  Livermore  National  Laboratory. 

Below  is  a  summary  of  CFDRC’s  Response  and  Comments  to  this  report  (in  the  frames  are 
extracted  relevant  fragments  from  the  original  report): 


Coventor  Designer 

Commercial  Mesher 

CFDRC  Micromesh 

speed  of  defining  geometry 

good 

adequate 

very  slow;  volumetric 
extrusion  is  the  bottleneck 

CFDRC  Response: 

-  The  “volumetric  extrusion”  is  relatively  slow  in  Micromesh  only  for  non-rectangular  shapes  in 
masks,  when  millions  of  fine  voxels  must  be  used  on  large  areas.  If  layout  masks  contain 
rectangular  shapes  only,  then  voxels  are  not  used  at  all,  and  the  3D  model  generation  is  much 
faster. 


meshing  speed 

good 

good 

slow.  However,  it  takes 
more  man-hours  to  "get 
ready  to  mesh"  than  it 
takes  to  numerically 

generate  the  mesh. 

-  A  totally  new  algorithm  and  procedures  for  prism-hex  mesh  generation  in  Micromesh  have 
been  implemented  and  completed  in  4*'^  qtr  2002.  The  meshing  speed  is  more  than  2  orders  of 
magnitude  faster  now!  Moreover,  the  resulting  mesh  quality  is  much  better,  where  the 
minimum  angle  in  all  triangles  can  be  controlled  (typically  now  the  min.  angle  is  >  10  deg). 


API  available 

no 

simple  enough  to  reverse 

no 

engineer 

API  available 


Is  it  possible  to  access  the  meshing  routines  from  library  with  a  well-documented 
Application  Program  Interface?  Failing  that,  is  it  possible  to  call  the  mesher  from 
the  command  line? 


-  Also  in  4*  qtr  2002,  CFDRC  completed  the  implementation  of  a  new  API  and  a  Python- 
based  scripting,  which  allows  access  to  all  the  Micromesh  functions,  including  importing 
layouts,  geometric  model  building,  and  all  the  meshing  routines.  The  documentation  is  being 
currently  prepared. 

Moreover,  it  is  possible  to  call  the  mesher  from  the  command  line,  with  parameters  like 
layout  file  (GDSII  or  CIF),  technology  file,  etc. 
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These  new  capabilities  allow  to  use  Micromesh  in  an  EDA-style  design  flow,  controlled  by 
external  tools. 


3.  Often,  it  has  been  necessary  to  manually  replace  GDSII  coordinates  with  "hand-typed"  coordinates, 

rounded  to  some  tolerance,  in  order  to  generate  an  acceptable  mesh.  We  do  not  know  if  the  problem 
has  been  with  the  GDSII  input  or  if  it  is  a  problem  with  Micromesh.  In  either  case,  it  would  be 
beneficial  to  have  a  "snap-to-grid"  capability  for  imported  geometries,  e.g.  round  all  geometry 
coordinates  to  a  specific  tolerance.  It  should  be  noted  that  we  recently  received  a  Micromesh 
upgrade  that  has  a  new  procedure  for  tolerancing,  this  has  not  been  evaluated. 

-  The  "snap-to-grid"  capability  for  imported  geometries,  that  is  rounding  all  geometry 
coordinates  to  a  specific  tolerance,  has  been  implemented  in  the  new  version  of  Micromesh 
(available  in  the  upgrade  package  sent  to  LLNL  on  1/20/03).  This  facilitates  fixing  of  the  “bad” 
GDSII  layouts  which  have  traces  that  should  be  joined,  but  are  not  quite,  etc. 


6.  The  "cut  and  paste"  capability  of  Micromesh  is  limited  in  that  it  only  works  for  a  single  layer,  i.e.  it 

is  difficult  to  replicate  a  multilayer  geometry. 


-  The  capability  to  “cut  and  paste”  or  “copy  and  paste”  for  multiple  layers  at  the  same  time 
has  been  implemented,  and  is  currently  available  in  Micromesh  (since  March  2003). 

New  licenses  of  CFDRC  software,  including  Micromesh,  were  provided  with  no  charge  to  Univ. 
of  Cincinnati  (Ranga  Vemuri),  Univ.  of  Illinois  (Andreas  Cangellaris),  University  of  Delaware 
(Fouad  Kiamilev’s  group),  and  University  of  Washington  (Richard  Shi’s  group)  for  NeoCAD 
applications  and  interoperability  studies. 

The  overall  progress  of  the  CFD-Micromesh  tool  during  the  DARPA  NeoCAD  program  is 
illustrated  below. 
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Examples  of  building  time  of  high-resolution  3D  model  and  mesh 
in  progressing  Micromesh  technologies 


Micromesh  1.4  (2000) 

Micromesh  2.0  (2001) 

Micromesh  2003 
(May  2003) 


32-iines  interconnect 
structure 

120  minutes 

(7200  sec.) 


2  minutes 
(120  sec.) 

2  seconds 


64-lines  interconnect 
structure 


00  (impossibi 

with  voxels) 

10  minutes 
(600  sec.) 

10  seconds 


Figure  62.  Overall  progress  of  the  CFD-Micromesh  tool  during  the  DARPA  NeoCAD  program. 


New  Micromesh  Applied  to  Another  NeoCAD  Benchmark 

The  newest  Micromesh  version,  including  all  the  newly  developed  technologies  mentioned 
above,  was  applied  to  another  NeoCAD  Benchmark:  UW-SOI-QVCO  placed  on  the  DaVinci  site 
by  University  of  Washington  (UW).  This  is  an  SOI  quadrature  voltage-controlled  oscillator 
(QVCO)  with  three  inductors  placed  at  varying  distances  to  test  magnetic  coupling  effects  on 
VCO  phase  noise  performance.  The  layout  proposed  by  UW  is  presented  below. 
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5.6GHz  Quadrature  LC  Oscillator  Layout 


inductors 


M22 


Surrounding 

inductor 


Figure  63.  NeoCAD  Benchmark  layout:  UW-SOI-QVCO  from  University  of  Washington. 


The  GDSII  file  provided  by  UW  in  the  benchmark  set  was  imported  into  Micromesh,  and  a  3D 
model  was  built,  which  is  shown  below. 


Figure  64.  The  UW-SOI-QVCO  benchmark  3D  model,  automatically  built  from  the  GDSII  layout, 

using  the  new  Micromesh  capabilities. 
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Next,  using  the  new  Micromesh,  we  generated  automatically  a  3D  structured  mesh,  fragment  of 
which  is  presented  below.  It  is  a  quite  regular,  hexagonal  mesh,  which  is  quite  convenient  for 


Figure  65.  A  3D  structured  mesh  generated  automatically  by  Micromesh  from  UW  layout. 


Summarizing,  Micromesh  became  during  the  NeoCAD  Program  a  very  useful  tool  for  automatic 
generation  of  complex  3D  models  and  computational  mesh  from  layout  formats  used  in  EDA 
community  (GDSII  and  GIF).  Micromesh  has  been  successfully  used  by  several  NeoCAD  teams 
and  has  become  an  indispensable  component  of  full-wave  EM  modeling  software  in  several 
NeoCAD  design  tools  (e.g.  Orora,  LLNL,  U.  Wash.,  U.  Pitt.,  SUNY,  and  others). 


6  VHDL-AMS  COMPACT  MODELS  OF  O/E  DEVICES 

Task  led  by  University  of  California,  Santa  Cruz,  PI:  Prof  Steve  Kang 

During  the  NeoCAD  OE-IDEAS  effort,  we  developed,  implemented,  and  tested  compact  models 
for  various  opto-electronic  (O/E)  devices,  with  the  aim  of  creation  of  a  model  library  for  mixed- 
level  mixed-technology  simulations.  The  new  models  were  demonstrated  through  example 
VHDL-AMS  simulations.  The  following  O/E  devices  were  addressed  and  tasks  accomplished. 

6.1  VCSEL  and  Transmitter  (VCSEL  and  a  CMOS  driving  circuit) 

-  Developed  a  behavioral  VHDL-AMS  model  of  VCSEL; 

-  Implemented  the  model  into  SMASH  from  Dolphin  Integration; 
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-  Tested  the  VCSEL  model  by  time  domain  simulation  of  different  devices; 

-  Simulated  a  simple  transmitter  consisting  of  CMOS  driving  circuit  and  VCSEL. 

Automation  of  0/E  behavioral  model  parameters  extraction: 

-  Finished  the  linear  part  of  VCSEL  behavioral  model  parameter  extraction  using  Python; 

-  Implemented  a  nonlinear  behavior  model  parameter  extraction  algorithms  for  Python. 

6.2  Metal-Semiconductor-Metal  (MSM)  Photodetector 

-  Developed  a  behavioral  VHDL-AMS  model  of  MSM  photodetector; 

-  Implemented  the  model  into  hAMster  provided  by  AnSoft; 

-  Validated  the  model  by  comparing  simulation  results  with  experimental  data. 

6.3  Raman  Amplifier 

-  Developed  a  behavioral  VHDL-AMS  model  of  Raman  Amplifier; 

-  Implemented  the  model  into  hAMster  from  AnSoft; 

-  Validated  the  model  by  comparing  simulation  results  with  experimental  data. 

6.4  Erbium  Doped  Fiber  Amplifier  tEDFA) 

-  Developed  a  behavioral  VHDL-AMS  model  of  EDFA; 

-  Implemented  the  model  into  hAMster  from  AnSoft; 

-  Validated  the  model  by  comparing  with  the  experimental  data. 

A  detailed  description  of  the  derived  compact  models,  with  references,  and  representative 
VHDL-AMS  simulation  results  are  presented  in  Appendix  A  at  the  end  of  this  report. 


7  DESIGN  AND  FABRICATION  OF  TEST  STRUCTURES  AND  VALIDATION  FOR 
RF  INTERCONNECTS  AND  OPTO-ELECTRONIC  SYSTEMS 

In  this  task.  University  of  Delaware  and  Cornell  University  groups  designed,  fabricated,  and 
measured  test  0/E  and  mixed-signal  VLSI  structures  to  validate  models  and  extract  necessary 
physical  parameters. 

7.1  Design,  Fabrication,  and  Testing  of  OE  and  RF  chips  and  PCBs 

This  task  was  primarily  performed  by  University  of  Delaware  (UD),  led  by  Dr.  Fouad  Kiamilev. 
UD  has  been  focusing  on  providing  experimental  circuits  and  measured  data  for  the  verification 
and  validation  of  CFDRC  new  modeling  and  simulation  capabilities  in  the  OE-IDEAS  project. 

First  year  accomylishments: 

•  Completed  the  design,  fabrication  and  testing  of  a  testing  chip  using  TSMC  0.25um 
CMOS  technology.  The  chip  contains  arrays  of  VCSEL  drivers,  optic  receivers,  LVDS 
drivers  and  receivers. 

•  Transferred  the  design  database  to  Micromesh  for  subsequent  processing. 
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Second  year  accomplishments: 

•  Generated  Mayo  lOG  system  MCM  design  database  in  DXF  format  for  CFD  to  transfer 
the  database  to  Mieromesh  and  to  simulate  the  MCM  electrical  interconnects 
performance. 

•  Performed  experimental  study  of  the  signal  integrity  performance  of  the  electrical  traces 
on  DARPA  RATS  main  board. 

•  Performed  the  theoretical  and  simulation  studies  of  Via-Stub  resonance  phenomenon  on 
system  PCBs. 

•  Designed  and  implemented  low  power  Current  Mode  Logic  electrical  interconnects  in 
Peregrine  0.5um  Silicon  on  Sapphire  process.  Testing  is  in  process. 

Graphics: 

OE-IDEAS  Test  Chip  fist-year  accomplishment): 


Figure  66.  Chip  with  Bonding  Wires. 
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LVDS Driver-Receiver  Testins: 


VDDA3.3V 


Digital  Control 
Bias  line  1 
Bias  line  2 

Digital  In  ' 
enable 


TX 


IDACI  IDACI 


VSSA  OV 
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Features: 

•  Independent  and  digitally  controlled  bias 
current  and  modulation  current  settings 

•  Power  down  c^ability 
Current  steering  for  reduced  SSN 


Layout  and  photogr^h 


LVDS  Driver-Receiver  Testing: 


Signal 

Source 
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Mayo  lOG  Opto-Electronic  System  Multi-Chip  Module  (MCM) 

find-year  accomylishment): 


Figure  67.  Mayo  lOG  system  MCM  layout,  in  Cadence-Allegro  (left)  and  CFD-Micromesh  (right). 


60 


Voltage  at  the  far  end  of  the  aggressor  line  (V)  Voltage  at  the  far  end  of  the  aggressor  line  (V) 


Measurement  of  Signal  Intesrity  Properties  on  DARPA  RATS  Board  Electrical  Traces 

(2’'  -year  accomplishment) 


_  „  White:  Victim  line  _  _ 

To  Oscilloscope  ^ - 1  | - ►To  Oscilloscope 


Gray:  Aggressor  line 
To  Signal  Source  M - 1  ~ 


►To  Oscilloscope 


art  of  DARPA  Rats  Board 
iShowing  Electrical  Traces 


2.0 

1.5- 

1.0 


i  I*  ^  ^  Ji  |N  . 


0.5- 

0.0 

-0.5- 

-1.0- 

-1.5- 


Time  (ns)  ^ooo 


-  Aggressor  Line 


0.0  ^ 


E 

-0.1  ."C 


Cross-talk  on  TRL3  with  near  end  of  the  victim  terminated  by  50n 


-0.5  ^ 

§ 


Time  (ns)  2000 


Time  (ns) 
-  Aggressor  Line 


'"1“  tK^T 


--0.5- 

u 

5-1.0- 


Cross  Talk  Measurements 


*  J.  J*|  *  *  iK 


V'JwV’wVVwwV 

Cross-talk  on  TRL3  with  near  end  of  the  victim  open 


Time  (ns)  2000 

Blue:  aggressor 
Red:  victim 


1000  Time  (ns)  2000 


1000  Time  (ns)  2000 


61 


Voltage  at  the  far  end  of  the  victim  line  (V)  Voltage  at  the  far  end  of  the  victim  line  (V) 


Dual-Rate  Current  Mode  Logic  Electrical  Interconnects  (2nd-vear  accomplishment) 


Dual  Rate  CML 
Electrical 
Interconnects  were 
implemented  in 
Peregrine  0.5um 
SOI  process; 
Testing  is  on-going 


Figure  68.  Dual-Rate  Current  Mode  Logic  Electrical  Interconnects  test  chip. 


7.2  Characterization  and  modeling  of  on-chip  interconnect  with  realistic  loads  and  drives 

Full  3D  Maxwell's  solver  cannot  be  directly  applied  to  the  complex  interconnect  problem  in 
mixed-technology  IC  due  to  its  huge  complexity.  Any  computational  algorithm  for  interconnect 
analysis  has  to  be  linear  or  0(n),  where  n  is  the  number  of  interconnect  lines  in  the  system.  The 
constraint  mainly  comes  from  the  increase  in  n  also  following  the  growth  in  Moore's  Law,  and 
hence  we  cannot  otherwise  anticipate  future  computing  machines  to  handle  the  complexity. 
Algorithms  of  0(n)  almost  dictates  that  the  analysis  will  be  based  on  lumped  or  spectrum 
analysis.  At  Cornell  University,  the  S-parameters  approach  with  generalization  to  large  signals 
has  been  attempted  for  interconnect  coupling.  Interconnect  structures  have  been  fabricated  in 
Cornell  Nanofab  with  the  initial  measurements  presented  in  the  figures  below. 


62 


Figure  69.  Interconnect  structures  fabricated  in  Cornell  Nanofab. 


Transmission 


S12  =  0.2  > 
S21  =  0.2  > 


S12  =  0.008  <  97.43 
S21  =  0.008  <  98.14 


Ginax=  -23.492  dB 


Reflection 


k  =  116.432 


Figure  70.  Results  of  initial  measurements  of  the  interconnect  structures 
fabricated  at  Cornell  University. 


Transient  Analysis  of  Nonlinear  Microwave  Circuits  Using  Small-Signal  S-Parameters 

The  number  of  carrier  frequencies  limits  the  application  of  Harmonic-Balance  (HB)  method, 
which  has  matured  as  the  workhorse  for  nonlinear  microwave  circuit  analyses.  A  newly 
proposed  frequency-domain  method  in  allows  more  accurate  description  of  the  frequency 
dependence  than  that  of  the  Root  model.  It  can  also  be  applied  to  strongly  nonlinear  circuit 
analysis,  where  the  conventional  frequency-domain  methods,  such  as  power  series  and  Volterra- 
series  analysis,  are  not  appropriate.  However,  the  total  number  of  frequencies  that  can  be 
considered  is  still  limited. 
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Microwave  circuits  intended  for  telecommunication  applications  are  expected  to  handle  one  or 
more  carriers  modulated  with  non-null  information  signals,  i.e.,  finite  bandwidth  excitations. 
The  bandwidth  issue  is  more  pronounced  in  CDMA  and  the  emerging  ultra-wide-band  (UWB) 
communication  systems.  Furthermore,  waveform  engineering  for  the  linearization  of  power 
amplifiers  in  CDMA  systems,  waveform  analysis  of  high-efficiency  power  amplifiers  (such  as 
class  F)  and  UWB  systems  require  accurate  transient  analysis  and  simulation  methods  for 
nonlinear  microwave  circuits,  which  was  considered  not  necessary  a  decade  ago.  The  broadband 
nature  of  those  problems  makes  the  available  methods  difficult  to  use. 

Here  we  propose  a  general  method  for  nonlinear  circuit  transient  analysis  based  on  measured 
small-signal  scattering  parameters.  The  method  does  not  suffer  from  the  limitations  of  the 
aforementioned  approaches.  The  procedure  directly  constructs  the  large  signal  response  from  the 
small-signal  scattering  parameters.  The  method  is  applied  to  the  simulation  of  a  high  electron 
mobility  transistor  (HEMT)  power  amplifier.  Reasonable  agreement  between  the  simulated  and 
measured  results  is  obtained.  This  method  does  not  have  limitations  on  the  number  of  input 
carrier  frequencies  or  the  total  number  of  frequencies.  It  is  expected  to  be  a  useful  and  efficient 
tool  in  waveform  engineering  applications. 


Figure  71.  Approximation  of  input  waveform. 
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Figure  72.  Measured  and  simulated  results,  (a)  Output  power.  The  square  is  the  simulated  points, 
and  the  triangles  are  measured  points;  (b)  Simulated  output  voltage  waveforms. 
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For  more  details  on  this  topie,  see  the  attached  paper  from  Cornell  University  “Transient 
Analysis  of  Nonlinear  Microwave  Circuits  Using  Small-Signal  Scattering  Parameters”,  by 
Pingshan  Wang,  Valeriy  S.  Kaper,  Shealy  J.  Richard,  and  Edwin  C.  Kan,  submitted  to  IEEE 
MTT-S  Conference,  June  2003. 

Interconnect  modeling  of  mixed-signal  integrated  circuits 


A  small  signal  s-parameter  transient  analysis  method  for  high-speed  interconnect  with  CMOS 
inverters  as  transceivers  is  developed.  The  input  signal  is  approximated  as  the  summation  of 
small  amplitude  step  functions.  The  linear-microwave-circuit  theory  is  used  in  the  simulation 
algorithm  with  the  measured  scattering  parameters  of  the  nonlinear  active  devices  and  the 
dispersive  interconnect  lines.  Thus,  no  frequency  dependent  information  is  lost.  The  self- 
consistent  method  may  also  be  used  to  analyze  the  small  signal  crosstalk  simultaneously. 
Simulation  results  agree  well  with  experimental  results. 


A  typical  interconnect  representation  is  shown  in  Figure  73.  The  interconnect  line  is  generally  a 
linear  system  (except  the  nonlinear  transmission  line  interconnect  proposed  previously  in  this 
NeoCAD  project  as  soliton  interconnect),  despite  some  operating  voltage  level  related 
phenomena.  Its  frequency  dependent  characteristics  are  faithfully  and  conveniently  described 
using  small-signal  s-parameters,  Sl,  shown  in  Figure  73.  The  driver  and  receiver  are  nonlinear 
devices.  Their  small  signal  s-parameters  depend  on  the  bias  conditions.  As  a  result,  their  s- 
parameters  across  the  full  bias  space.  So  and  Sr,  will  be  used  in  the  following  simulation  and 
analysis. 


z 


L 


Figure  73.  An  interconnect  circuit  with  an  input  power  source,  an  input  network  with  s-parameter 
Sin,  and  an  output  network  with  S„ut.  Zs  is  the  source  impedance.  The  input  and  output  networks 
provide  connections  between  the  power  source  and  the  load.  Tj  and  El  denote  the  source  and  load 

reflection  coefficient  in  the  frequency  domain. 

A  transient  analysis  needs  to  find  the  circuit  response  to  an  arbitrary  input  signal,  shown  in 
Figure  74.  We  approximate  the  input  signal  as  the  summation  of  small  amplitude  step- 
functions,  each  with  a  different  time  delay,  shown  in  Figure  74. 
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The  ith  step  function  is 

(1) 

i=l 

where  n  is  the  number  of  discrete  steps.  The  Fourier  transform  of  the  ith  step-function  is, 

Vs,,  )  =  (^^.,,+1  -  ){7rd (o))  +  — )  exp(  -jat, )  (2) 

JO) 

The  response  of  the  circuit  to  the  i*  small  step  input  pulse  can  be  found  using  the  microwave 
circuit  theory  through  the  signal  flow  graph  in  Figure  75  provided  that  So  and  Sr  are  specified 
using  the  algorithm  in  Figure  76. 
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Figure  75.  (a)  A  signal  flow  graph  of  the  circuit  in  Fig.  1.  A  Zo  (reference  impedance  of  the  s- 
parameters)  is  splitted  from  the  source  impedance,  Z*,  in  Fig.l.  The  rest  of  the  sources  impedance 
is  combined  with  the  input  network,  and  represented  by  A,„.  The  separation  is  not  necessary,  but 
convenient  for  analysis  formula  simpliflcations.  (b)  Simplified  signal  flow  graph  by  cascading  the  s- 

parameters  of  all  the  elements  in  (a). 


In  Figure  75,  we  have 


Therefore,  we  have 


With  the  ineident  and  reflected  waves  given  in  equation  (5),  the  voltages  and  currents  at  any  port 
in  circuit  Figure  73  can  be  obtained  easily  with  the  help  of  the  scattering  parameters  of  the 
relevant  sub-network.  An  inverse  Fourier  transformation  gives  the  voltages  and  currents  at  the 
interested  time  moment. 
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Step  approximation  of 
input  waveform 


Calculate  the  circuit  response  to  the  small  step  input 

pulse 


Figure  76.  A  flow  chart  of  the  frequency-domain  algorithm 

A  new  bias  point  for  the  next  small  step  input  calculation  at  the  following  time  moment  is  then 
obtained  by  the  summation  of  the  circuit  responses  to  all  the  previous  small  step  inputs  at  the 
interested  time  moment.  Using  the  algorithm  in  Figure  76,  the  transient  response  of  the  circuit 
can  then  be  obtained. 

The  arbitrary  input  signal  in  Figure  74  can  also  be  approximated  as  the  summation  of  narrow 
square  functions  as, 

U(0  =  XU,(«C (6) 

i  =  l 

where  n  is  the  number  of  discrete  steps.  The  Fourier  transform  of  the  ith  square  signal  is, 

K,  ((o)  =  K,,-  —  (- exp(-7^yt,+i )  +  exp(-7^yt,. ))  (7) 

JO) 

The  system  response  to  this  square  signal  depends  on  the  bias  conditions  of  the  nonlinear  devices 
determined  by  their  responses  to  all  the  previous  square  signals  at  the  moment  when  this  signal 
appears.  If  the  time  step  (the  width  of  the  square),  U+i-U,  is  not  large,  the  amplitude  of  each 
frequency  component  in  (7)  is  small.  The  circuit  operation  condition  will  not  change  much  as 
well.  Therefore,  the  system  is  an  instantaneous  linear  system  to  this  square  pulse  though  pulse 
amplitude  seems  large.  Similar  equations  and  algorithm  can  be  used  for  the  transient  analysis. 
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In  the  concluding  phase  of  the  interconnect  modeling,  we  solidified  our  programs  for  large-signal 
reconstruction  from  S-parameter  measurements.  Benchmarks  of  typical  data  paths  containing 
nonlinear  amplification  elements  are  compared  with  full  Maxwell  solver,  reduced  solver,  and 
SPICE-level  large-signal  simulation.  Also,  we  have  consolidated  our  effort  in  using  permalloy 
and  undergrid  decoration  of  transmission  lines  (XL)  and  coplanar  waveguides  (CPW)  to  achieve 
broadband  interconnect  with  small  dispersion  in  MHz-20GHz  ranges.  This  technology  not  only 
makes  the  design  of  complex  interconnect  network  more  tractable  from  decreased  cross  talk  and 
more  manageable  large-signal  behavior,  but  also  gives  the  available  interconnect  simulation 
algorithm  better  handles  to  reduce  the  complexity.  If  no  localization  assumption  can  be  applied 
through  either  the  fdament  or  partial  inductance  concept,  the  complex  geometry  of  full-chip 
interconnect  (more  than  10^  geometry  objects)  usually  renders  the  modeling  work  intractable. 
This  is  one  problem  that  we  cannot  hope  to  be  solved  by  better  computing  capabilities,  because 
to  our  present  knowledge  better  computers  can  only  be  made  from  more  complex  interconnect 
network,  and  the  simulation  has  to  be  linearly  scaled. 


Permalloy  Loaded  Transmission  Line  for  High-speed  Interconnect  Applications 

Self  inductance  and  mutual  inductance  of  global  interconnect  in  deep  sub-micrometer  CMOS 
technologies  are  crucial  factors  of  integrated  circuits  performance.  Self  inductance  depends  on 
the  current  return  path,  which  is  unknown  prior  to  extracting  and  simulating  the  entire  nonlinear 
circuits.  A  large  number  of  the  neighboring  interconnects  need  to  be  considered  for  accurate 
estimation  of  the  mutual  inductance.  Consequently,  extraction  of  inductance  is  the  present 
bottleneck  in  high-speed  interconnect  modeling  and  design. 

The  partial  inductance  method  is  the  general  practice  for  such  problems,  evaluating  inductance 
for  each  wire  segment  by  an  imaginary  loop  with  infinity.  Furthermore,  building  a  complete 
inductance  matrix  for  all  interconnects  on  an  integrated  circuit  (IC)  is  impractical  because  this 
procedure  will  result  in  a  large,  tightly  coupled  network  representation.  Unlike  capacitance 
matrices  that  can  be  truncated  to  represent  only  localized  couplings,  simply  discarding  remote 
mutual  inductances  can  produce  an  unstable  equivalent  circuit  mode.  Based  on  this  increased 
complexity  and  what  is  known  from  full-chip  capacitance  extraction,  any  attempt  at  full-chip 
inductance  extraction  requires  some  form  of  localization  via  windowing,  which  is  hardly 
applicable  to  long  wires.  Techniques  used  in  capacitance  extraction,  such  as  library  construction 
and  analytical  formulas,  cannot  be  directly  applied  to  inductance  extraction.  Several  other 
methods,  such  as  return-limit  loop  inductance,  the  2x  mutual  inductance  screening  rule,  and  the 
susceptance  matrix,  have  been  proposed.  These  methods  have  their  own  specific  limitations 
related  to  the  penetrating  nature  of  magnetic  field  coupling.  Therefore,  the  ultimate  solution  to 
the  interconnect  inductance  extraction  and  modeling  is  to  eliminate  the  mutual  magnetic  field 
coupling  and  the  dependence  on  the  current  return  path  for  the  self  inductance.  This  can  result  in 
a  diagonalized  inductance  matrix  with  elements  accurately  determined  by  the  partial  inductance 
method.  We  show  that  the  incorporation  of  patterned  permalloy  in  interconnects  provides  such  a 
path.  Introduction  of  high  frequency  magnetic  material  increases  the  self  inductance. 
Nonetheless,  it  is  demonstrated  that  inductance  can  have  positive  effects  on  the  performance  of 
VLSI,  which  suggests  that  increasing  inductance  in  a  controllable  way  can  be  beneficial  for 
overall  system  considerations. 
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The  figure  below  shows  the  test  structures  fabricated  using  a  CMOS  compatible  process.  Lines 
with  or  without  patterned  permalloy  are  designed  for  direct  comparison  on  self  and  mutual 
inductance.  Lines  with  patterned  permalloy  bars  are  2pm  wide,  7pm  long  and  around  200nm 
thick.  A  patterned  Tantalum  (Ta)  film  of  ~20  nm  thickness  serves  as  an  adhesion  and  migration- 
barrier  layer.  The  distance  between  two  adjacent  permalloy  bars  is  2pm.  Aluminum 
transmission  lines  are  5pm  wide.  The  coupling  section  of  the  two  lines  that  are  parallel  to  each 
other  is  3,705pm  long.  The  separation  between  the  parallel  transmission  lines  is  10pm.  These 
dimensions  are  all  layout  dimensions  subject  to  lithography  variations,  except  for  the  thickness. 
Identical  comparison  structures  without  patterned  permalloy  were  also  fabricated.  The  ground 
aluminum  plate  is  ~8pm  below  the  top  metals  separated  by  Si02  deposited  at  ~460nm/minute 
using  plasma-enhanced  chemical  vapor  deposition  (PECVD).  The  porosity  in  PECVD  oxide  is 
high  enough  to  relax  the  residual  stress  that  no  serious  wafer  warping  is  observed. 


Figure  77.  (a)  Schematic  of  the  magnetic  coupling  test  structures;  (b)  part  of  the  fabricated  test 

structures 


The  comparison  structures  are  designed  to  have  a  strong  magnetic  coupling,  but  small  capacitive 
coupling.  The  estimated  coupling  capacitance  and  mutual  inductance  for  the  comparison 
structure  are  \%.9fF  and  367p//,  respectively.  The  structures  with  permalloy  have  slightly  larger 
coupling  capacitance  and  essentially  the  same  mutual  inductance. 

After  performing  measuring  pad  de-embedding  and  equivalent  circuit  elements  extraction 
procedures,  the  figure  below  shows  the  measured  S12  while  port  3  and  4  are  left  open. 
Asymptotically  similar  curves  were  obtained  if  port  3  and  port  4  were  connected  to  a  50  Ohm 
load  through  several  sub-miniature  adapter  (SMA)  elbow  connectors.  The  measured  for  the 
lines  with  patterned  permalloy  structures  is  essentially  from  the  noise  background.  The 
comparison  structure  shows  very  small  coupling  at  low  frequencies.  The  coupling  grows  along 
with  the  increase  of  frequency  as  expected  for  the  test  structures.  It  is  clear  that  patterned 
permalloy  has  effectively  reduced  the  magnetic  coupling  for  more  than  60  dB,  and  the  mutual 
inductance  is  virtually  eliminated. 
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For  more  details  on  this  topic,  see  the  paper  from  Cornell  University  “Permalloy  Loaded 
Transmission  Line  for  High-speed  Interconnect  Applications”,  by  Pingshan  Wang  and  Edwin 
Kan,  submitted  to  IEEE  MTT-S  Conference,  June  2003. 

Patent  application  has  been  filed  by  the  Cornell  University  group. 


8  NEOCAD  TEST  STRUCTURES  AND  BENCHMARKS  SIMULATIONS 
8.1  NeoCAD  Benchmarks  from  the  DaVinci  Netbook  website 

In  May  2002,  CFDRC  downloaded  all  the  materials  from  the  DaVinci  Netbook  website 
containing  the  benchmark  cases  suggested  for  NeoCAD  participants.  We  built  3D  models  and 
ran  3D  EM  simulations  of  several  selected  benchmark  cases.  This  effort  is  reported  below. 

As  a  first  test  structure  we  took  HRL-Mayo  InP  test  chip  “ic_hrl_tcl”,  posted  on  the  DaVinci 
web  site.  On  that  chip  Mayo  designed  a  variety  of  test  structures,  including  coupled  Line  Pairs 
(LP),  Coplanar  lines  (CP),  and  a  Microstrip  line  (MS)  structure  (see  Figure  79).  The  substrate 
was  700-micron  thick  semi-insulating  InP.  Two  levels  of  gold  interconnect  were  separated  by  a 
polyimide  dielectric  layer  1.4  microns  thick.  Metal- 1  lies  on  the  substrate  and  is  0.55  microns 
thick.  Metal-2  is  on  top  of  the  polyimide  dielectric  and  is  2  microns  thick. 
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Figure  19.  Layouts  of  several  test  structures  on  the  Mayo  test  chip  “ic  hrl  tcl”. 


The  icjirljcl  GDSII  layout  file  was  imported  into  CFD-Micromesh  and  analyzed  with  respect 
of  the  most  appropriate  test  structures  for  the  initial  simulations  with  NeoCAD  tools  at  CFDRC. 
As  the  first  test  case,  we  selected  the  CPI  structure  (Figure  79),  that  is  Coplanar  line  of  Metal  2, 
with  Spacing  Sum  and  Width  Sum.  The  CPI  layout,  imported  in  CFD-Micromesh,  is  shown  in 
Figure  80,  together  with  the  computational  mesh.. 


After  setting  up  in  CFD-Micromesh  the  technological  data  (thicknesses  and  positions  of  all  the 
layers)  of  the  ic_hrl_tclch\^,  a  S-D  computational  model  was  generated  automatically  by  CFD- 
Micromesh.  An  example  of  structured  grid  is  shown  in  Figure  80  below. 
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Figure  80.  Layout  of  the  CPI  structure,  imported  in  CFD-Micromesh  from  GDSII  file,  and  the 
computational  grid  generated  automatically  by  CFD-Micromesh  from  the  layout. 


Next,  we  performed  a  series  of  static-electric  and  transient  full-wave  electromagnetic  simulations 
with  ACE+,  to  extract  the  basic  parameters  (resistance,  I-V  characteristics,  S-parameters)  of  the 
structure,  to  be  compared  with  the  measured  data  posted  on  the  DaVinci  site,  in  the  Mayo-Bl 
directory. 

Before  performing  a  full- wave  transient  electromagnetic  simulation  of  the  CPI  structure,  a 
steady-state  electric  simulation  of  the  line  was  performed,  to  calculate  the  resistances  of  all  the 
lines.  An  example  snapshot  of  these  results  is  shown  below  (with  3.3  V  applied  to  the  ground 
pads  A  and  C).  Both  the  line  resistance  values  agree  very  well  with  the  values  measure  at  Mayo 
-  see  the  table  below.. 
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Figure  81.  An  example  snapshot  of  steady-state  electric  simulation  results  of  CFD-ACE+,  with  3.3 
V  applied  to  the  ground  pads  A  and  C,  and  zero  volts  to  all  the  other  pads. 


Line  Resistances: 


Calcuiated 
with  CFD-ACE+ 

Measured 
at  Mayo 

B-E  Line 

19.76  Q 

19.86  Q 

A-D  Line  (ground) 

5.12  Q 

5.15  Q 

Next,  a  full-wave  transient  electromagnetic  simulation  of  a  signal  pulse  propagation  on  the  CPI 
serpentine  coplanar  waveguide  was  performed  using  CFDRC’s  EMAG  simulator.  A  Gaussian 
pulse  with  amplitude  3.3  V  was  applied  between  the  input  (pad  B)  and  the  ground  pads  (A  and 
C).  Several  selected  snapshots  of  that  transient  simulation  are  shown  in  the  figures  below,  as  well 
as  resulting  time-dependent  voltage  and  current  waveforms,  and  calculated  from  them  S- 
parameters  converted  to  frequency  domain. 
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Figure  82.  Simulated  Gaussian  Pulse  Propagation  in  Time  in  the  CPI  Benchmark  Structure. 
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Figure  83.  Voltage  and  Current  Monitors  at  the  CPI  line  beginning  (red)  and  end  (green). 


Figure  84.  Measured  and  Calculated  S-Parameters  for  the  CPI  Test  Structure. 
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20G  Waveguide  Benchmark  from  NCSU 


Another  structure  taken  from  the  DaVinci  web  site  for  CFDRC  software  validation  was  the  20G 
structure  from  NCSU-B2B  directory.  The  structure  is  illustrated  in  Figure  85  below. 
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Figure  85.  The  20G  test  structure  from  NCSU-B2B:  basic  parameters  and  cross-section  photo. 

For  our  first  calculation  we  selected  the  structure  of  the  width  Wm  =  0.5  pm,  and  the  “short” 
length  Ls  =  400  pm.  A  structured  computational  mesh  built  automatically  in  CFD-Micromesh  is 
shown  in  Figure  86.  Using  this  mesh,  we  performed  full-wave  electromagnetic  simulations  with 
CFD-Maxwell. 
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Figure  86.  Computational  mesh  built  in  CFD-Micromesh  for  20G  structure  from  NCSU. 


The  computed  S-parameters  are  shown  in  Figure  87.  For  comparison,  the  S-parameters  measured 
by  NCSU  are  plotted  in  Figure  88.  Some  issues  may  be  observed.  First,  note  that  the  measured 
data  cover  much  lower  range  of  frequencies  (45  MHz  -  20  GHz),  while  the  calculated  range  of 
frequencies  is  9GHz  -  200  GHz.  This  is  due  to  the  fact  that  our  current  CFD-Maxwell  algorithms 
do  not  allow  to  calculate  long  simulated  time  intervals  due  to  very  small  computation  time  step 
(limited  by  CFL  number,  or  Courant  factor).  Therefore,  the  short  simulated  time  intervals 
translate  into  very  high  frequency  range  after  FFT  transformation.  This  difficulty  should  be 
overcome  after  implementation  of  the  ADI  algorithm,  reported  earlier  in  this  report. 

Second,  even  if  both  in  the  measurements  and  calculations  Sll  is  very  close  to  S22,  and  S12  is 
virtually  equal  to  S2 1 ,  the  values  obtained  from  experiments  and  calculations  are  very  different. 
In  the  measured  data  the  reflection  coefficients  (Sll,  S22)  are  high  (0.6... 0.9),  and  the 
transmission  coefficients  (SI 2,  S21)  are  much  lower  (0.2  to  0.3  for  higher  frequencies).  In  our 
calculated  results  it  is  the  opposite:  in  the  similar  frequency  range  (up  to  20  GHz,  or  2E+10)  the 
reflection  coefficients  (Sll,  S22)  are  lower  (less  than  0.5),  and  the  transmission  coefficients 
(S12,  S21)  are  much  higher  (more  than  0.9).  This  discrepancy  was  discussed,  also  with  the 
authors  of  the  measurements  from  NCSU  (Prof.  M.  Steer’s  group).  The  reflection  and 
transmission  depend  very  much  on  the  line  termination  impedance  (matching),  and  we  were  not 
sure  if  the  experimental  equipment  (microprobes  and  pads)  setup  had  not  affected  the  measured 
characteristics  (if  the  line  impedance  was  not  ideally  matched  on  both  ends,  the  reflection 
coefficients  Sll  and  S22  would  be  much  higher). 
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Figure  87.  S-parameters  computed  by  CFD-Maxwell  for  the  20G  structure. 
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Figure  88.  S-parameters  of  the  20G  structure,  measured  by  NCSU. 
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8.2  Simulation  of  Critical  Path  from  the  Mayo  “lOG”  System  MCM  Board 


In  the  first  year  of  NeoCAD,  after  extensive  discussions  with  participation  of  CFDRC,  U- 
Delaware,  U-Pittsburgh,  and  Mayo,  it  was  decided  that  the  'lOG'  FSOI  system  funded  under 
DARPA  VLSI-Photonics  Program  would  be  used  a  benchmark  problem  for  the  teams  working 
on  integrated  electronic/photonic  system  design.  A  sketch  of  the  lOG  System  is  depicted  in 
Figure  5.1. 

According  to  the  NeoCAD  goal  of  using  the  lOG  System  as  a  demonstrator,  we  decided  to  use 
the  MCM  (multi-chip  module)  board  of  the  lOG  System  as  an  example  for  our  electromagnetic 
simulations.  The  most  interesting  would  be  analysis  of  a  “critical  path”  on  the  MCM  board, 
which  may  limit  the  overall  performance.  Such  a  critical  path  has  been  identified  through 
consulting  with  Mayo,  UPitt,  and  UDelaware.  It  is  illustrated  graphically  in  the  following 
figures. 


Figure  89.  Sketch  of  the  “lOG”  FSOI  System,  accepted  as  the  NeoCAD  benchmark. 

The  MCM  board  is  at  the  bottom. 
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Figure  91.  Single  chip  layout  on  the  lOG  MCM  board. 


Blue  line  :  signal  l,  Red  line  :  signal_2 
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Figure  92.  Some  signal  traces  in  the  layer  “signal  !  chip_U4”.  The  C4out  signal  path  (two 
differential  lines)  has  been  selected  as  the  “critical  path”  for  first  electromagnetic  simulations. 


The  C4  “critical  signal  path”  from  the  above  layout  of  the  lOG  System  MCM  board  was  selected 
to  be  simulated  by  CFDRC’s  full- wave  electromagnetic  simulator,  to  extract  crucial  electrical 
parameters  (e.g.  S-parameters,  cross-talk,  delays). 

The  layout  of  the  full  C4  Path,  after  importing  into  Micromesh,  is  shown  in  the  figure  below. 
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Figure  93.  The  C4out  signal-path  layout,  and  a  3D  model  built  in  Micromesh. 
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Figure  94.  Mesh  details  of  the  C4out  Line,  with  static  differential  electric  potential  results. 


Before  performing  a  full-wave  transient  electromagnetic  simulation  on  the  critical  path  C4out,  a 
steady-state  electric  simulation  of  the  line  was  performed,  to  calculate  the  line  resistance.  A 
snapshot  of  these  results  is  shown  below  (with  -I-3/-3  V  applied  to  the  differential  lines  C4out  and 
N_C4out).  Both  the  line  resistance  values  calculated  with  ACE-i-  were  less  than  3%  different 
from  the  values  extracted  by  Cadence  Allegro  tool. 
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Figure  95.  Snapshot  of  the  steady-state  electric  simulation  results  from  CFD-ACE+,  with  +3  V  and 
-3  V  applied  to  the  differential  lines  C4out  and  N_C4out. 


Next,  a  full-wave  transient  electromagnetic  simulation  of  a  signal  pulse  propagation  on  the 
critical  path  C4out  was  performed  using  the  CFDRC’s  EM  AG  simulator.  A  differential  signal 
(Gaussian  pulse)  of  +0.4  and  -0.4  volt  was  applied  to  the  input  of  the  differential  lines  C4out  and 
N_C4out.  Several  selected  snapshots  of  that  transient  simulation  are  shown  in  the  figures  below. 
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Figure  96.  Several  selected  snapshots  of  that  transient  simulation  for  a  differential  signal  (Gaussian 
pulse)  of  +0.4  and  -0.4  volt  applied  to  the  input  of  the  differential  lines  C4out  and  N_C4out. 

The  Impedance  and  Propagation  Delay  calculated  for  the  C4out  line  by  CFDRC  full-wave 
EMAG  solver  were  within  20%  difference  range  as  compared  to  the  values  extracted  by  Cadence 
Allegro  tools. 


9  CONCLUSIONS,  COMMERCIALIZATION,  TECHNOLOGY  TRANSFER 

9.1  Conclusions  and  Summary  of  Achievements 

This  report  has  presented  the  current  status  and  results  of  the  OE-IDEAS  project  within  DARPA 
NeoCAD  Program,  aimed  at  developing,  demonstration,  and  validation  of  new  modeling 
concepts  and  tools  in  the  following  domains: 

•  new,  efficient  techniques  of  extraction  of  parasitics  as  well  as  new  methods  of  full-wave 
EM  modeling  of  parasitics  in  electrical  interconnects,  where  high-frequency 
electromagnetic  effects  do  not  allow  for  using  simplified  SPICE-type  models  anymore; 

•  new,  comprehensive  CAD  tools  for  integrated  opto-electronic  (0/E)  circuits  and  systems, 
which  would  allow  designers  to  analyze  and  optimize  the  full-path  electrical-optical  signal 
characteristics,  at  all  the  stages  of  the  optically  communicating  circuits. 

The  major  achievements  of  the  program  have  already  been  outlined  in  the  Introduction  section, 
and  here  are  presented  again  in  brief: 
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•  First  in  the  U.S.  commercial  software  for  high-fidelity  multi-mode  simulation  of 
Lasers  (VCSELs,  edge-emit.)  and  Photodetectors,  O’SEMI  from  CFDRC: 

-  new  Optical  Gain  Model,  with  quantum-physics  and  thermal  phenomena; 

-  implemented  Effective  Index  Method  for  Edge-Emitting  Laser  optics; 

-  new  Tunneling  Current  Model,  based  on  quantum  model  of  heterojunction; 

-  currently  installed  and  used  by  AFRL/WL,  ARE  (Adelphi,  MD),  NRL,  UCSB, 

UC  Santa  Cruz,  U.  Florida,  U.  Arizona,  U.  Waterloo,  Vanderbilt  Univ.,  and  others; 

-  inquiries  from:  Sandia  National  Labs,  LENT,  McMaster  University,  et  al. 

•  New  capabilities  in  CFDRC’s  Full-Wave  Electromagnetic  Simulators 
for  RF  and  Optical  Interconnects: 

-  Signal  monitoring  and  postprocessing:  V(t),  I(t),  Z(w),  S-param(w); 

-  Electromagnetic  filament  model  implemented  and  initially  tested; 

-  CFDRC  full- wave  EM  solvers  offered  to  NeoCAD  Community,  and 
used  by  Univ.  of  Delaware,  U.  of  Washington,  and  U.  of  Cincinnati. 

.  CFDRC’s  Full-Wave  3D  Simulator  EMAG  coupled  with  SPICE. 

•  Demonstrated  Photonics  Applications  and  provided  Tutorials  for  EMAG  users: 

1)  Optical  waveguide  coupling; 

2)  Optical  ring  resonator; 

3)  Photonic  crystal  waveguide; 

4)  Photonic  bandgap  2D  matrix. 

•  Simulated  several  NeoCAD  benchmark  structures  and  demonstrated  results. 

•  CFD-Micromesh  -  a  unique,  automated  3D  mesher  from  IC  layouts  -  enhanced,  offered  to 
NeoCAD  Community;  new  tutorials  provided. 

-  Micromesh  became  during  the  NeoCAD  Program  a  very  useful  tool  for  automatic 
generation  of  complex  3D  models  and  computational  mesh  from  layout  formats  used  in 
EDA  community  (GDSII  and  CIF). 

Micromesh  has  been  successfully  used  by  several  NeoCAD  teams  and  has  become  an 
indispensable  component  of  full-wave  EM  modeling  software  in  several  NeoCAD 
design  tools  (e.g.  Orora,  LENT,  U.  Wash.,  U.  Pitt.,  SUNY,  and  others). 

•  At  the  University  of  Delaware: 

-  Designed,  fabricated,  tested  and  measured  an  optical  transceiver  chip,  and  built  a 
3D  model  (with  CFD-Micromesh). 

-  Characterized  and  built  a  3D  model  (with  CFD-Micromesh)  of  the  Mayo  lOG  system 
multi-chip  module  (MCM). 

-  Formulated  the  transmission-line  model  and  the  7i-model  of  via-stub  resonance  on 

system  PCB’s;  with  CFDRC  simulated  the  model  performance  . 

-  Designed,  fabricated,  and  tested  Dual-Rate  Current  Mode  Logic  electrical  10 ’s. 

•  At  Cornell  University: 

-  Established  the  quantum  transport  model  hierarchy  in  density  gradient,  quantum 
hydrodynamic,  full  Wigner,  and  NEGF  on  large  barrier  tunneling  problems.  The 
benchmark  problems  are  RTD  and  MOSFET  gate  tunneling. 
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-  Established  large-signal  transient  analysis  from  S-parameters  on  nonlinear  circuits. 
The  benchmark  problems  include  high-power/high  frequency  transistors,  nonlinear 
transmission  lines,  and  buffer-insertion  interconnect. 

-  Proposed  and  analyzed  pulse-based  interconnect  architecture  and  circuit 
implementation,  based  on  the  large-signal  transient  analysis. 

-  Designed  and  verified  the  use  of  patterned  permalloy  to  control  the  on-chip  magnetic 
fields  for  accurate  partial  inductance  extraction. 

-  Designed  and  verified  the  use  of  undergrids  to  control  CMOS  interconnect 
dispersion  on  transmission  lines  and  co-planar  waveguides. 

•  At  the  University  of  California  Santa  Cruz: 

Developed,  implemented,  and  tested  compact  models  for  mixed-level  mixed- 

technology  simulation,  and  demonstrated  in  VHDL-AMS  for  the  following  OE 

devices: 

-  VCSEL  and  transmitter  (including  driver  CMOS  circuit) 

-  Metal-Semiconductor-Metal  (MSM)  photodetector: 

-  Raman  Amplifier: 

-  Erbium  Doped  Fiber  Amplifier  (EDFA): 

-  Automation  of  0/E  behavioral  model  parameters  extraction  from  CFDRC  3D  high- 
fidelity  3D  simulator,  for  VCSEL  and  MSM  photodetector. 


9.2  Commercialization  and  Technology  Transfer 

Commercialization:  This  project  resulted  in  several  components  with  strong  commercial  value, 

including: 

•  Advanced  3D  full-wave  electromagnetics  solver,  EMAG,  solving  full  Maxwell  equations, 
has  been  developed  and  implemented  as  a  integrated  module  in  CFD-ACE+  Multiphysics 
software.  It  is  currently  marketed  and  distributed  by  ESI-Group,  but  CFDRC  has  the  rights 
for  adaptation  and  applications  for  DoD  and  US  defense  industry  projects. 

•  Micromesh  -  An  automated  geometric  modeler  and  mesher,  customized  for  electronic 
circuits  applications.  CFDRC  will  commercialize  it  directly  as  a  separate  product,  and  as  a 
module  in  3'^'^  party  EDA  software  (Orora,  Cadence, ...). 

•  O’SEMI  -  CFDRC-owned  software  tools  for  modeling  and  design  of  semiconductor 
electronic  and  optoelectronic  devices.  CFDRC  is  working  with  selected  academic  teams  and 
DoD  labs  on  software  validation  and  qualification  for  commercial  use.  O’SEMI  will  be 
marketed  to  US  industry  in  the  form  of  customized  tools,  projects,  and  ultimately  as  a 
licensable  software  tool. 

•  CFD-Maxwell  -  CFDRC-owned  flexible  FDTD  solver  of  full  3D  Maxwell  equations,  with 
powerful  capability  of  body-conforming  and  unstructured  mesh  solutions,  for  high-fidelity 
electromagnetic  and  optical/photonic  design  problems.  It  is  being  evaluated  by  academic 
teams,  and  will  be  used  for  projects  with  DoD  labs,  defense  industry,  for  commercial 
industry  projects,  and  a  customized  version  as  a  licensable  technology. 
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Technoloev  Transfer:  The  Technology  Transfer  resulting  from  the  OE-IDEAS  project  includes: 

•  25  publications  (by  all  team  members) 

•  1  patent  (from  Cornell  University) 

•  CFDRC  NeoCAD  3D  EM  Solver  (and  Micromesh)  to  be  interfaced  with  Analog/RF 
Design  Flow  from  Orora  Design  Technologies 

•  New  proposed  SBIR  project:  “CAD  Framework  for  Coupled  EM+  Circuit  Simulation” 
has  been  accepted  (in  March  2004)  for  funding  by  DARPA/MTO. 

In  the  new  “Coupled  EM+Circuh”  SBIR  project,  CFDRC  and  Dr.  Vikram  Jandhyala  from 
University  of  Washington  will  be  working  closely  with  Orora  Design  Technologies  towards 
bringing  the  new  CAD  framework  and  tools  to  RF/analog  circuit  designers  and  EDA  market. 
The  existing  Orora  tools  (Arsyn,  Arcap  and  Arana)  can  be  used  standalone,  or  used  with  tools 
from  Cadence,  Mentor  Graphics,  and  other  third-party  tools,  to  form  a  complete  analog/RF 
design  flow  (see  Figure  97  below). 
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Figure  97.  Analog/RF  Design  Flow  with  related  Orora  tools,  and  the  new  Coupled  EM+Circuit 
tools  to  be  implemented  in  the  next  CFDRC  SBIR  project  (starting  Spring  2004). 


CFD-Micromesh:  Fast  3D  Model  and  Mesh  Builder  from  Layouts 

The  new  CFD-Micromesh  software,  greatly  enhanced  during  this  project,  for  automatic  3D  solid 
model  generation  and  meshing  from  electronic  IC  layouts,  has  a  very  big  and  promising 
commercial  potential.  The  popular  layout  description  formats,  CIF  and  GDSII,  are  imported 
easily  into  CFD-Micromesh,  allowing  coupling  the  software  with  several  commercial  EDA 
design  tools.  The  3D  model  building  and  meshing  in  CFD-Micromesh  is  very  fast  and  controlled 
by  a  user-friendly  graphical  user  interface  (GUI),  as  well  as  by  external  Python  scripts  or 

command  lines.  The  automatically  generated  three-dimensional  device  models,  using  the 
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structured  or  unstructured  3D  computational  mesh  in  DTP  format,  can  be  imported  directly  into 
CFDRC's  simulator  CFD-ACE+  for  multi-physics  simulations,  or  to  other  NeoCAD  EM  solvers, 
after  translation  to  the  neutral  mesh  format  also  developed  under  the  NeoCAD  Program. 

The  Cornell  University  team  (headed  by  Prof  Edwin  Kan)  will  try  to  secure  more  external 
funding  to  continue  the  interconnect  technology  development  and  modeling. 

•  Both  non-linear  transmission  line  (NLTL)  and  permalloy  technology  have  found 
applications  in  the  DARPA  Chip-scale  Atomic  Clock  (CSAC)  project  under  PM  of  A1 
Pisano. 

•  The  large-signal  reconstruction  from  small-signal  responses  is  continuing  development  with 
more  test  cases  and  the  AWE  method  (asymmetric  waveform  estimation). 

•  The  pulse-wave  interconnect  has  published  its  first  architectural  paper,  and  the  future 
silicon  implementation  is  uncertain  at  this  stage.  The  work  related  to  pulse-wave 
interconnect  needs  more  evidence  before  commercialization  can  be  considered. 

•  The  permalloy  transmission  line  work  is  under  patent  application. 

The  University  of  California  Santa  Cruz  team  (headed  by  Prof.  Steve  Kang)  plans  to  share  the 
developed  compact  models  of  OE  devices  as  broadly  as  possible. 

•  As  a  first  step,  we  plan  to  apply  the  models  for  a  study  of  CD  writing  using  blue  laser 
diodes. 

•  For  system  studies  and  potential  commercial  applications,  Prof  Kang’s  PhD  student  (Shu 
Wu)  will  work  with  Prof  Ken  Pedrotti  of  UCSC  to  perform  feasibility  studies  sponsored  by 
National  Semiconductor  and  the  State  of  California. 


Interaction  of  CFDRC  Team  with  DoD  Laboratories 


Since  the  beginning  of  the  NeoCAD  program,  CFDRC  has  established  collaborative  partnerships 

with  three  DoD  Laboratories  in  areas  directly  related  to  this  project: 

•  AFRL  Wright  Patterson  AFB.  Contact:  Dr.  Thomas  Nelson,  Tactical  Lasers  Branch, 
Sensors  Directorate.  CFDRC  has  delivered  ACE-i-  software  to  AFRL  for  modeling 
optoelectronic  devices  (VCSELs,  Photodetectors).  AFRL  arranged  a  DoD  wide  training  of 
ACE-i-  at  AFRL  WE  in  Dayton,  OH,  on  Oct  1-3,  2002.  Seven  people  from  Air  Force  and 
Army  attended  the  workshop.  One  of  the  key  areas  of  interest  was  electronic,  thermal,  and 
optical  modeling  of  VCSELs. 

•  ARL  Adelphi  MD.  Contact  Dr.  Brian  Riely,  Optoelectronics  Branch.  CFDRC 
representative  visited  ARL  in  June  2002,  and  Dr.  Riely  attended  ACE-i-  training  course  at 
CFDRC  in  Huntsville  in  Aug  2002.  Dr.  Riely  is  already  using  ACE+  at  ARL.  The  major 
interest  of  this  collaboration  is  modeling  of  optoelectronic  devices  (VCSELs)  and  effects  of 
optical  emission  on  quantum  wells. 

•  NRL  Washington  DC:  Dr.  Felix  Buot,  Power  Electronics  Branch.  Dr.  Buot  spent  one  year 
sabbatical  at  CFDRC  working  on  enhancements  of  quantum  effects  models  in  ACE-i-  SEMI 
module,  and  on  modeling  power  semiconductor  devices.  CFDRC  has  installed  ACE-i-  at  NRL 
and  Dr.  Buot  is  using  the  code  for  Navy  research  projects. 
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Interaction  with  US  Industry 


During  the  course  of  this  project  CFDRC  has  established  technical  interaction  contacts  with 
several  US  industry  partners  including: 

•  Rockwell  Science  Center,  Thousand  Oaks  (Dr.  Rong  Chow)  joint  project  for  analysis  of 
wide  bandgap  semiconductor  devices,  and  packaging  problems  of  RF  power  electronics. 

•  Mixed  Technology  Associates  (MTA),  Palo  Alto,  CA,  (Dr.  Daniel  Yergeau):  adaptation  of 
ACE+  modeling  tools  for  MTA  TCAD  device  simulation  tools  for  RF  devices. 

•  LSI  Logic  (Ken  Doniger):  joint  analysis  of  a  test  set  consisting  of  three  problems: 

1)  A  spiral  inductor  with  substrate:  calculate  inductance  and  Q  of  the  device. 

2)  Three  parallel  conductors  of  trapezoidal  cross-section:  calculate  the  capacitance  as  a 

function  of  the  trapezoidal  angle. 

3)  A  thin  metal  serpentine  structure:  extract  resistance  and  capacitance 

•  IBM,  T.  J.  Watson  Research  Center,  Yorktown  Heights,  NY,  (Dr.  Steven  Koester),  a 
member  of  DARPA  TEAM  Program:  joint  submission  of  a  new  SBIR  proposal 
“Computational  Analysis  of  Electromagnetic-Noise  Tolerance  of  Nanoelectronic  Devices” 
to  DARPA. 

•  Eastman  Kodak  Company,  Integrated  Materials  and  Microstructures  Lab,  Rochester,  NY, 
(Dr.  Edward  P.  Furlani,  Senior  Principal  Scientist):  since  beginning  of  March  2003,  CFDRC 
has  been  providing  software  and  support  for  simulation  of  photonic  structures  (waveguides, 
photonic  crystals),  and  for  optoelectronic  semiconductor  device  simulation  (VCSELs, 
photodetectors).  Kodak  plans  to  use  CFDRC’s  software  to  design  and  optimize  new 
optoelectronic  devices. 

•  Agilent  R&D  Center,  San  Jose,  CA,  (Dr.  Darek  Burak):  presentation  of  CFDRC’s 
Optoelectronic  software  and  discussion  of  future  collaboration  and  software  applications. 

•  Lumileds  Lighting  Advanced  Labs,  San  Jose,  CA,  interested  in  modeling  multiple  quantum 
well  light  emitting  diodes  and  design  of  anti-reflection  surface  features.  Got  O’ SEMI  trial 
license,  software  sales  is  being  negotiated. 

Several  companies  (e.g.  Samsung  from  Korea,  NEC  from  Japan)  have  expressed  interest  in 
using  CFD-ACE-i-  electromagnetic  tools  (developed  under  this  project)  for  modeling  full- wave 
electromagnetic  and  parasitic  effects  in  VLSI  analog/RF  ICs  and  interconnects.  CFDRC  has 
developed  example  demonstration  cases  and  tutorials,  which  should  lead  to  software  sales  and/or 
additional  R&D  projects. 

Interaction  with  Academia 


More  than  600  academic  users  have  been  using  ACE-i-  Multiphysics  software  worldwide.  The 
ACE-i-  EMAG  module  for  VLSI/RF  Integrated  Circuits,  however,  has  not  been  formally  released 
yet.  At  present  CFDRC  provides  a  "research"  version  to  few  selected  academic  partners.  Some  of 
them  have  committed  to  participate  in  joint  development  and  enhancements  of  ACE-i-  EMAG. 
Academic  teams  using  or  evaluating  the  CFDRC  electromagnetic  codes  include:  Cornell  Univ. 
(Prof  E.  Kan),  SUNY  (Prof  A.  Pacelli),  Georgia  Tech  (Prof  G.  Cassinovi),  Univ.  Cincinnati 
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(Prof.  R.  Vemuri),  Univ.  of  Illinois  Urbana-Champaign  (Prof  A.  Cangellaris),  Univ.  Washington 
(Dr.  P.  Nikitin),  and  Univ.  of  Delaware  (Dr.  F.  Kiamilev). 


The  newly  developed  O’ SEMI  package  is  currently  used  by  several  academic  groups,  including: 
two  licenses  at  UC  Santa  Barbara  (Prof  Larry  Coldren,  Prof  Joachim  Piprek),  UC  Santa  Cruz 
(Prof  S.  Kang),  University  of  Arizona,  ECE  Department/Optical  Sciences  (Prof  M.  Neifeld), 
Univ.  of  Florida  (Prof  Fan  Ren’s  group,  collaborating  with  ARE),  Wilfrid  Laurier  University  in 
Waterloo,  Ontario  (Prof.  M.  Wartak),  Univ.  of  Texas  in  Arlington  (Dr.  B.  Borovic),  Vanderbilt 
University  (Prof.  Ron  Schrimpf,  Prof.  Greg  Walker). 


Talks  /  Presentations  (not  listed  in  the  published  References) 

•  DARPA/MTO  NeoCAD  Kickoff  Meeting  in  Malibu,  California  (Aug.  29-30,  2001): 
OE-IDEAS  Project  presentation  was  given  by  the  PI,  Marek  Turowski. 

•  DARPA/MTO  NeoCAD  Program  Review  Meeting  in  San  Antonio,  TX  (Feb.  27-28,  2002): 
OE-IDEAS  Project  presentation  was  given  by  the  PI,  Marek  Turowski. 

•  NeoCAD  OE-IDEAS  Project  Review  by  AFRL  Monitors:  Norman  Bernstein,  Paul  Cook, 
(Rome  Labs),  and  Thomas  Nelson  (WP-AFB),  at  CFDRC  in  Huntsville,  AL  (May  2,  2002): 
presentations  were  given  by  Marek  Turowski  and  Andrzej  Przekwas. 

•  DARPA/MTO  TEAM/NeoCAD  2002  Summer  Review  in  Providence,  RI  (Aug.  27-29, 
2002):  OE-IDEAS  Project  presentation  was  given  jointly  by  the  PI,  Marek  Turowski,  and  a 
subcontractor,  Edwin  Kan  of  Cornell  Univ. 

•  Presentation  at  Agilent  R&D  Center  (San  Jose,  CA),  entitled  “Computational  Multiphysics 
Technology  at  CFDRC:  Applications  to  Optoelectronics  and  Photonics”,  by  Zbigniew 
Sikorski  (CFDRC),  on  Jan  29,  2003. 

•  Presentation  “Multiphysics  modeling  of  intracavity  VCSELs“  at  the  Photonic  West  2003 
Conference,  Session  Physics  and  Simulation  of  Optoelectronic  Devices  XI,  given  by 
Zbigniew  Sikorski  (CFDRC),  on  Jan  29,  2003. 

•  DARPA/MTO  NeoCAD  2003  Winter  Review  in  Arlington,  VA  (Feb.  19-20,  2003): 
OE-IDEAS  Project  presentation  was  given  by  the  PI,  Marek  Turowski. 

•  DARPA  TEAM/NeoCAD  2003  Summer  Review  in  Monterey,  CA  (Sep.  15-17,  2003): 
OE-IDEAS  Project  presentation  was  given  by  the  PI,  Marek  Turowski. 

Patents 

A  new  nonlinear  interconnect  scheme  for  integrated  circuits,  named  soliton  interconnect,  has 
been  conceived  at  the  Cornell  University  (PI:  Prof  Edwin  Kan),  using  microwave  circuit 
techniques  to  enhance  information  transaction  on  the  interconnect  network.  The  new  technology 
has  been  filed  for  patent  consideration. 
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9.3  Recommendations  and  Plans  for  Further  Research 


Based  on  the  experience  gained  in  this  project  we  have  identified  several  areas  for  further 
research,  model  development,  and  CFD-ACE+  software  improvement.  They  can  be  grouped  in 
following  categories: 

a)  Interfaces  with  EDA  Tools 

In  the  present  project  we  have  greatly  enhanced  CFD-Micromesh,  to  interface  microelectronic 
and  3D  EM  computational  domains.  The  code  is  rapidly  gaining  acceptance  and  interest  from 
broad  range  of  applications.  It  is  one  of  the  biggest  achievements  of  this  project.  Micromesh  has 
however  still  several  needs  for  system-level  design,  for  example: 

•  semiautomatic  block  structured  grid  generation, 

•  mixed  structured-unstructured  grids, 

•  adaptive  octree  grids, 

•  parametric  modeling  in  Micromesh, 

•  automated  construction  of  filaments  in  Micromesh, 

•  tighter  link  with  process  fabrication  and  layout  tools, 

•  boolean  operations  on  masks  for  multilayer  microstructures. 

b)  Automation  and  Parameterization  of  Model  Setup 

At  present,  CFD-ACE-l-  is  suitable  for  single  runs  of  a  specific  design  or  pre-arranged  multi-runs 
for  single  variation  in  e.g.  boundary  conditions.  There  is  a  strong  need  for  parametric  modeling, 
scripting,  optimization,  and  inverse  problem  solutions.  A  first  step  has  been  recently  made  at 
CFDRC  by  developing  a  Python  scripting  capability  for  our  tools.  We  hope  to  expand  this 
capability  to  more  applications.  Advanced  academic  and  industrial  users  should  be  able  to  link 
the  functions  and  modules  in  ACE+  modeling  dataflow.  Coupling  of  ACE+  with  an  optimization 
module  would  answer  requests  from  several  customers  who  have  requested  device  optimization 
capability. 
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This  Appendix  A  to  OE-IDEAS  Final  Report  is  eomposed  of  four  ehapters. 

In  chapter  I  we  present  VHDL-AMS  implementation  of  the  rate  equation  based  Vertical-Cavity 
Surface-Emitting  Laser  (VCSEL)  model,  discussion  of  the  time  domain  simulation  results  for 
different  kinds  of  VCSEL,  a  prototype  of  a  transmitter  comprised  of  a  VCSEL  and  a  CMOS 
driving  circuit,  and  the  model  parameter  extraction  issue. 

In  chapter  II  we  present  VHDL-AMS  implementation  of  a  Metal-Semiconductor-Metal 
photodetector  (MSM-PD)  model,  and  simulation  results. 

In  chapter  III  we  present  VHDL-AMS  implementation  of  an  Erbium  Doped  Fiber  Amplifier 
(EDFA)  model,  and  simulation  results. 

In  chapter  IV  we  present  VHDL-AMS  implementation  of  a  Raman  Amplifier  (RA)  model,  and 
simulation  results. 
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CHAPTER  I 


MODELING  AND  SIMULATION  OF  VCSEL  AND  TRANSMITTER 


1.1.  Model  description 

The  model  in  [Mena,  1999]  starts  from  spatially  dependent  rate  equations.  To  reduce  the 
computational  complexity  the  carrier  distribution  was  modeled  by  Bessel  series  and  mode  profile 
was  modeled  by  a  Laguerre-Gaussian  function  [Mena,  1999],  Though  integration  we  obtained 
spatially  independent  rate  equations  as  the  circuit-level  VCSEL  model,  which  is  a  series  of 
ODEs  (Ordinary  Differential  Equations).  The  complete  mathmatical  model  and  physical 
explanation  can  be  found  in  [Mena,  1999]. 

The  VHDL-AMS  implementation  of  this  VCSEL  model  is  straightforward.  After  the  port  map  of 
the  device  model  was  determined,  the  ODEs  can  be  directly  imbedded  into  the  model.  It  is 
should  be  noted  that  even  though  our  implementation  and  simulation  are  in  SMASH™  provided 
by  Dolphin  Integration,  after  some  modification  the  VHDL-AMS  VCSEL  model  can  be  used  by 
any  simulation  tool  supporting  VHDL-AMS.  According  to  the  ability  of  the  simulation  tool  we 
can  change  the  time  unit  of  the  model  parameters  to  help  achieve  the  convergence  in  simulation. 
Following  code  expresses  the  implementation  of  VCSEL  by  VHDL-AMS. 

library  IEEE; 

use  IEEE. math_real. all; 

entity  VCSEL  is 
generic! 

-parameters  of  VCSEL— 
wq:  real:=  1.60219e-19;  — electron  charge 
wgO:  real:=  3.0e4; 

port(  terminal  p,  n,  phout:  electrical); 
end  entity  VCSEI; 

architecture  Singlemode  of  VCSEL  is 
-quantities  and  variable  declaration 
quantity  lin  through  p  ;  -input  current 

-  use  break  to  set  the  initial  condition 
break  NO  =>  0.0; 

-set  variable  value 
gain:=... 

-  rate  equations 

NO'dot  ==  wita*lin/(wq)-NO/(wtaon)-...; 

-  output  laser  power 
pw:=  wkfO*  SO; 

end  architecture  singlemode; 
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1.2.  Simulation  results  and  discussion 


1.2.1.  VCSELs 

As  a  typical  demonstrative  simulation  example,  we  simulated  the  response  of  a  VCSEL  under 
two  types  of  input,  namely  square  wave  and  sinusoidal  wave.  Simulation  results  are  shown  in  Fig 
1-1. 


(a)  Square  wave  input 


(b)  Sinusoidal  wave  input 

Fig  1-1.  VCSEL  responses  under  different  inputs 


From  the  simulation  results  we  observe  that  the  typical  large  signal  response  of  VCSEL  is 
oscillation  followed  by  a  steady  state  output  which  follows  the  input  signal  wave  form.  The 
differences  lie  in  the  oscillation  amplitude  and  time.  Simulation  results  showed  the  oscillation  of 
laser  under  square  wave  input  was  much  severer  than  the  sinusoidal  input.  It  is  because  the  rise 
and  falling  edges  of  square  wave  are  much  sharper  than  the  sinusoidal  wave  the  suddenly 
inputted  carriers  need  more  time  to  be  distributed.  The  laser  diode  would  need  more  time  to 
reach  its  steady  state. 

Due  to  the  lack  of  the  experimental  data  for  large  signal  response  of  VCSEL,  we  chose  an  eye 
diagram  reported  in  [Gustavsson,2002]  to  verify  the  multimode  simulation  ability  of  our  model. 
This  850-nm  oxide-confinement  top-emitting  VCSEL  consists  of  an  n-GaAs  substrate,  a  22  pair 
Alo.giGao.ogAs-Alo.iiGao.gsAs  p-DBR  on  top,  a  33  pair  Alo.giGao.osAs-Alo.nGao.gsAs  n-DBR  on 
the  bottom,  three  6-nm  GaAs  quantum  wells  sandwiched  by  2  8-nm  Alo.3Gao.7As  confinement 
layer,  and  a  40-nm  Alo.98Gao.02As  layer  to  create  index-guiding  oxide  aperture.  This  kind  of 
VCSEL  is  preferred  light  source  for  short  distance  high  speed  fiber  communication  system 
[Man,  1998]. 

The  simulation  results  show  that  when  the  drive  current  is  well  above  the  threshold,  an  optical 
bump  occurred  in  the  turn-off  transient.  Due  to  the  strong  index-guiding  along  the  transverse 
direction  in  oxide-confined  VCSELs,  the  transverse  optical  mode  is  confined  in  a  small  oxide 
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aperture.  As  a  result  the  carriers  in  the  central  region  are  depleted  much  faster  than  other  regions 
due  to  stronger  central  stimulated  recombination  processes.  When  the  injection  current  source  is 
suddenly  turned  off  the  carriers  in  the  peripheral  region  diffuse  to  the  central  region,  which  could 
make  the  VCSEL  turn  on  again.  This  spatial  hole  burning  (SHE)  effect  results  in  interaction 
between  the  optical  mode  and  carrier  profdes.  As  a  result  the  VCSELs  have  a  poor  turn  off 
characteristic.  The  simulation  results  in  Fig  1-2.  show  good  agreement  with  experimental  data  in 
[Gustavsson,  2002]  (Fig  1-10.  and  Fig  1-11). 
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(b)  drive  current  is  just  above  threshold 


Figure  1-2.  Simulation  results  vs.  experimental  data 


1.2.2  Transmitter 

As  we  mentioned  earlier  VHDL-AMS  is  a  good  choice  to  do  mixed-mode  simulation.  We  used  a 
simple  transmitter  to  show  this  advantage.  Even  though  the  simulation  was  not  done  for  a  real 
transmitter,  we  believe  that  after  reasonable  modification  of  the  parameters  we  can  use  this 
process  to  simulate  real  systems. 

The  schematic  of  a  driving  circuit  is  borrowed  from  [Chen,  1997].  The  VCSEL  is  a  typical  index- 
guided  InGaAs  VCSEL.  In  the  simulation  process  0.35  pm  technology  file  was  used.  Fig  1-3  are 
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the  simulation  results  of  driving  cireuit  and  transmitter.  In  Fig  1-3  (a),  V(INR)  represents  the 
input  signal  and  V(OUT)  is  the  output  of  the  driving  eircuit.  Under  a  200  MHz  1  V  to  -3.5  V 
input  voltage  signal  the  driving  eireuit  generated  a  200  MHz  2  mA  to  8  mA  eurrent  signal.  For 
speed  eonsideration  the  threshold  eurrent  of  the  VCSEL  is  ehosen  to  be  about  1.8  mA  to  make  it 
operated  above  threshold.  When  2  mA  eurrent  is  inputted,  the  VCSEL  would  generate  a  weak 
laser  signal  whieh  eould  not  be  deteeted  by  the  reeeiver.  The  peak  output  power  is  about  2  mw 
and  the  difference  between  the  peak  and  bottom  is  around  1.8  mW  which  is  not  desirable  in 
optical  communication  system.  By  introducing  a  larger  gain  VCSEL  such  as  thin-oxide-aperture 
we  can  get  a  better  transmission.  The  whole  simulation  procedure  suggests  that  because  of  its 
various  compatibility  VHDL-AMS  is  a  very  good  choice  to  do  the  optoelectronic  system 
modeling  and  simulation. 
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Fig  1-3.  Simulation  results  of  the  transmitter 
1.3.  Model  parameters  extraction 
1.3.1  Basic  flow  chart 


During  the  last  two  decades  much  research  has  been  done  on  modeling  and  simulation  research 
on  optoelectronic  devices  and  systems.  The  circuit-level  model  is  more  efficient  in  system  level 
simulation  compared  to  device-level  model,  which  requires  great  details  of  device  physics. 
Usually  this  kind  of  model  needs  multidimensional  analysis  of  spatial  behavior,  which  is  not 
suitable  for  system  simulation.  Another  advantage  of  circuit-level  model  is  that  it  can  be 
implemented  into  conventional  electronic  simulation  tool  such  as  SPICE.  An  important  before 
the  model  implementation  is  model  parameter  extraction  from  experimental  data  or  high  fidelity 
device-level  simulation. 


For  automatic  model  parameters  extraction  for  circuit-level  VCSEL  model,  we  chose  the  high 
fidelity  simulation  tool  CFD-ACE+.  By  taking  the  advantage  of  the  compatibility  of  Python  with 
CFD-ACE+,  we  can  take  Python  as  the  platform  to  implement  our  algorithm.  By  importing 
additional  numerical  enhancement  module,  Numpy,  the  implementation  can  be  done  more  easily 
than  directly  implementing  in  C.  The  main  procedure  is  shown  below. 
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Step  1.  Set  up  device-level  VCSEL  model  in  CFD-ACE-i-. 

Step  2.  Get  L-I  and  I-V  curve  in  different  temperature  through  Simulation  Manager. 

Step  3.  Extract  out  the  parameters  of  circuit-level  VCSEL  model. 

Step  4.  Put  the  parameters  extracted  into  circuit-level  VCSEL  model. 

Step  5.  Simulate  the  VCSEL  by  VHDL-AMS  simulation  tool. 

Step  6.  Compare  the  simulation  results  with  the  output  of  CFD-ACE-i-. 

Step  7.  Repeat  step  2  to  step  until  the  simulation  results  from  different  simulation  tools  are 
close  enough. 

The  corresponding  flow  chart  was  shown  in  Fig.  1-4. 


Fig  1-4.  Flow  chart  of  automatic  parameters  extraction  for  circuit-level  VCSEL  model 

procedure. 
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1.3.2  Linear  part  extraction 


1.3.2.1  Extraction  methodology 

In  order  to  extract  model  parameters  from  experimental  data  or  high  fidelity  simulation  tool  we 
need  to  use  curve  fitting  and  optimization  algorithms.  At  the  beginning  we  can  choose  least 
square  method  for  linear  fitting.  We  start  from  the  function  of  output  light  power  Pq 

Po=7(n(/-4-/<,//(n)  (1-1) 

where  T  is  temperature,  I  is  injected  current,  1*  is  threshold  current,  p  is  the  differential  slope 
efficiency  and  Lff  is  the  offset  current. 

The  offset  current  could  be  an  empirical  function  such  as  polynomial  function. 

V(0  =  E<i,r'  (1-2) 

1=0 

The  temperature  is  related  to  the  total  input  power,  IV,  and  output  light  power,  Pq,  as  following 
format. 


T  =  T,+{IV-P,)R„-t,,—  (1-3) 

at 

where  Rth  is  the  VCSEL’s  thermal  impedance,  Xth  is  a  thermal  time  constant,  To  is  the  ambient 
temperature,  and  V  is  the  laser  voltage.  Under  dc  conditions,  the  dT/dt  term  disappears.  Thus, 
from  the  resulting  equation  it  is  clear  that  (IV-  Po)  models  the  power  dissipated  in  the  VCSEL. 
Here  we  assume  that  any  power  not  carried  in  the  optical  output  is  dissipated  as  heat  in  the 
device.  From  L-I  and  I-V  curve  we  found  output  light  power  is  relatively  small  compared  with 
the  IV  term  so  we  can  neglect  Po  in  Eq.  (1-3).  The  simplified  relation  between  T  and  IV  is  given 
below. 


T  =  T,+IVR,,  (1-4) 

If  we  put  (1-2)  and  (1-4)  into  (1-1)  and  after  some  combination,  we  can  find  (1-1)  now  becomes 
the  following. 


Po=riiT) 


n 


'^{a,r+b,r) 


(1-5) 


Obviously  if  we  choose  n  as  a  large  number  it  would  be  more  accurate  but  also  more 
complicated.  Here  we  choose  n=4  for  both  accuracy  and  efficiency.  So  far  the  VCSEL  model  in 
CFD-ACE-i-  can  not  simulate  the  roll  over  of  L-I  characteristic  of  VCSEL,  and  we  can  set  ai  -  a4 
to  zero.  Fig  2-5.  shows  the  simulation  results. 
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Fig  1-5.  The  simulation  results  using  extracted  parameter  and  data  from  CFD-ACE. 
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CHAPTER  II 


MODELING  AND  SIMUALTION  OF  MSM  PHOTODETECTOR 


2.1.  Model  description  and  simulation  results 

2.1.1.  Numerical  model 

2.1.1.1. Model  setup 


The  behavior  of  photo-generated  carriers  in  MSM  photodetector  can  be  described  by  Poisson’s 
equation  and  carrier  continuity  equations  that  are  given  below 


V  =  -q{P  -n  +  N^-N^) 

q 


q  at 

Jn  =qPnnE„+qD„Vn 

Jp  =  qPpPEp-qDpVp 


(2-1) 


where  Rn/p  and  Gn/p  are  the  electron/hole  recombination  and  generation  rate,  \|/  is  the  electrostatic 
potential,  other  parameters  have  their  usual  meaning.  The  optical  generation  rate  is  assumed  to 
be  gaussian  for  the  impulse  response  and  is  given  as 


abs 


►  exp 


(  t  ^2^ 

(x-Xq) 

2 

exp 

2 

) 

1  ) 

(2-2) 


where  Wp  is  the  peak  optical  power  density  incident  on  the  surface,  hv  is  photon  energy,  aabs  is 
the  absorption  coefficient.  The  carrier  recombination  rate  can  be  modeled  using  Shockley-Hall- 
Read  (SHR)  theory.  After  assuming  that  excess  carrier  density  exceeds  the  intrinsic  carrier 
concentration  we  can  get  a  simple  form,  rJlxn  and  p/2Tp. 

We  want  to  reduce  the  complexity  for  circuit-level  model  extraction.  Since  most  of  MSM 
photodetectors  are  working  under  a  bias  voltage  which  is  larger  than  the  flat-band  voltage  we  can 
use  En=Ep=\[//d  instead  of  Poission’s  equation.  Above  equations  plus  infinite  recombination, 
(n,p(0,t)=n,p(L,t)=0),  boundary  conditions  complete  our  mathematical  model  for  MSM 
photodetector. 
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2.1.1.2.Simulation  results 


Fig.  2-1.  Photocurrent  in  MSM  photodetector 


The  model  above  can  be  solved  by  using  the  finite  difference  method.  We  chose  a  GaAs  based 
MSM  photodetector  to  describe  the  typical  response  of  a  MSM  photodetector  under  a  light 
impulse.  The  simulation  results  are  shown  below. 


Fig  2-1  shown  the  photocurrent  generated  by  the  device.  From  the  figures  we  can  see  that  since 
the  electrons  are  much  faster  than  holes  the  tail  of  the  holes  photocurrent  is  mostly  due  to  the 
holes. 


2.1.2.  Circuit-level  model 

2.1.2.1. Model  setup 

2.1.2.1.1.  Steady  state 

In  dc  steady  state  under  constant  illumination,  the  photocurrent  is  proportional  to  the  electron- 
hole  pair  generation  rate  as  shown  below 


108 


(2-3) 


where  Iph  is  the  photocurrent,  w  is  the  area  ratio  (area  transparent  to  light/  area  exposed  to  light), 
a  is  the  absorption  coefficient,  d  is  the  thickness  of  the  absorption  layer,  r\  is  the  quantum 
efficiency,  hv  is  the  photon  energy,  p  is  the  incident  light  power 

As  the  bias  voltage  is  raised  the  behavior  of  MSM  photodetector  can  be  treated  for  two  cases 
wherein  the  bias  voltage  is  lower  or  higher  than  the  flat-band  voltage.  When  the  bias  voltage  is 
lower  than  the  flat-band  voltage  the  potential  barrier  of  the  holes  is  lowered  by  the  increased  bias 
voltage,  which  results  in  rapid  increase  of  the  current,  when  the  bias  voltage  reached  the  flat- 
band  voltage,  the  hole  barrier  approaches  its  limiting  magnitude,  Schottky  barrier  height.  Further 
increase  in  voltage  will  cause  the  current  to  increase  slowly.  Thus  we  use  the  following  two 
formulas  to  model  the  photocurrent. 


^ dc  ^0 


^  dc  ^ ph 


qVdVps-V)  A 

e  -I  +l^j,-Ae'',V<Vf 

V  J 

r  A) 

l  +  Ae  ''  \,V>V, 


FB 


V 


(2-4) 

(2-5) 


FB 


J 


where  Idc  is  the  total  current,  V  is  the  bias  voltage,  Vfb  is  the  flat-band  voltage,  T  is  the 
temperature,  K  is  the  Boltzmann  constant,  n.  A,  B,  lo  are  fitting  parameters. 


We  implemented  this  model  using  VHDL-AMS  as  following. 


-  The  DC  model  of  MSM_PD  definition  begins 


LIBRARY  DISCIPLINES; 

LIBRARY  IEEE; 

USE  DISCIPLINES.ELECTROMAGNETIC_SYSTEM.ALL; 
USE  IEEE.MATH_REAL.ALL; 


-entity  declaration. 

ENTITY  MSMDC  IS 
GENERIC  (  T  :  REAL:=300.0;-k 
d  :  REAL:=1.0e-4;-cm 


): 

PORT(TERMINAL  p,m:  ELECTRICAL; 

quantity  ph  :real);  -Interface  ports. 
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END; 


--architecture  declaration. 

ARCHITECTURE  behav  OF  MSMDC  IS 
-quantity  declarations, 
constant  q  :real:=1 .6e-19;-C 
constant  he  :real:=6.6262e-34;-Js 
constant  k  :real:=1.3807e-23;-J/k 

QUANTITY  Iph  :  REAL; 

QUANTITY  lO  :  REAL; 

QUANTITY  vjn  ACROSS  i_out  THROUGH  p  TO  m; 

BEGIN 

ph==0.9823e-6/1 0.0**1 .5;-w 

lph==w*(1 .0-R)*(1 .0-exp(-1 .0*alfa*d))*eta*q*ph/hc/miu; 
IO==lph/(exp(q*Vfb**2/n/k/T)-1.0); 

if  v_in<Vfb  use 

i_out==IO*(exp(q*v_in*(2.0*Vfb-v_in)/(n*k*T))-1 .0)+lph*A*exp(-1.0*B/v_in); 
else 

i_out==lph*(1 .0+A*exp(-1 .0*B/v_in)); 
end  use; 

END; 

2.I.2.I.2.  Transient  response 

The  impulse  response  of  MSM  photodetector  can  be  modeled  as  a  finite  Fourier  sine  series 


^  2  .  f 

Jp'Y  0 

„=1  n;r  \  2  J 

(2-6) 

(2-7) 

where  Jp  is  fitting  parameter,  t  is  transit  time. 


The  responses  of  other  inputs  can  be  described  as  the  convolution  of  Jpuise  with  the  incident  light 
wave  function. 


J  other  =  Jpuise  *  PW)  =  J  p  ‘Y. - M  Jp  Sn  “  ^0  )PKto 

n=l 


(2-8) 


where  ph(t)  is  incident  light  wave  function. 


no 


2.1.2.2.Simulation  results 


To  demonstrate  the  model  we  did  simulation  corresponding  to  the  experiment  data  in  reference 
[Krishnamurthy,  2000]  and  [Aliberti,  2002],  Fig  2-2  corresponded  to  the  Fig  4  in 
[Krishnamurthy,  2000]  and  Fig  2-3  corresponded  to  Fig  1  in  [Aliberti,  2002],  The  parameters 
used  for  simulation  are  shown  in  the  figures.  Our  simulation  results  matched  those  in 
[Krishnamurthy.  20001.  [Aliberti.  20021. 


Fig  2-2.  Simulation  results  of  impulse  response  of  MSM-PD  in 
[Krishnamurthy,  2000] 
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Fig  2-3.  Simulation  results  of  impulse  response  of  MSM-PD  in  [Aliberti,  2002] 


It  should  be  noted  that  depending  on  the  ability  of  the  simulator  we  can  change  the  time  unit  of 
the  parameter  to  help  the  simulator  converge  to  a  correct  numerical  result. 

The  typical  cases  are  the  step  and  sinusoid  responses.  The  simulation  results  are  shown  in  Fig  2- 
4  and  Fig  2-5.  Because  there  is  no  published  experimental  data  available,  the  parameters  chosen 
in  simulation  do  not  correspond  to  the  real  devices.  The  transit  time  was  chosen  to  be  19.6  ps  in 
simulation  for  demonstration  purpose. 


Time  (s) 

Fig  2-4.  Simulation  result  of  the  MSM-PD  response  under  step 

function  input 
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CHAPTER  III 


MODELING  AND  SIMUALTION  OF  EDFA 


3.1.  Model  description 

Most  simply  we  can  model  EDFA  as  a  three  level  laser  system  as  Fig.  3-1.  So  we  can  write  rate 
equations  as  following  [Giles,  1991  ]: 


dN, 

dt 

dN, 

dt 

dN, 

dt 


=  -R,3N,  +  R3,N3  - 

=  W„N,-W,,N,-A,,N,+A3,N3 

=  Ri3Ni-R3iN3-A32N3 


(3-1) 


We  now  consider  the  steady  state  regime  where  the  populations  are  time  invariant.  From 
equations  (3-1)  we  can  get  steady  state  solutions 


N,=p 


ah 


b(a-i“  Rj3)  +  a\Vj2  +  Ri3''^32 


N,=p 


^13-^32  3.Wj2 


b(a  -h  Rj3  )  -h  a\Vj2  -t  Rj, A,, 


(3-2) 


where  p=Ni-i-N2+N3;  a==R3i-i-A32  and  b=W2i+A2i. 


When  a  signal  light  beam  with  intensity  L  (power  per  area)  at  wavelength  Xs  traverses  a  slice  of 
laser  medium  of  infinitesimal  thickness  dz  and  atomic  population  density  Ni  (lower  level)  and  N2 
(upper  level),  the  intensity  change  df  is  given  by 

dI.={<7i,(2.)Nj-<7,i(2.)N,  )l.dz  (3-3) 

where  012(^5)  and  a2i(Xs)  are  the  absorption  and  emission  cross  sections  of  the  laser  transition  at 
A-s,  respectively.  So  gain  g 

g  =  d<^.W)(l  +  D)-<J.(A)(l-D),  (3-4) 

where  D=(N2-Ni)/p=(2N2-p)/p,  with  D=-l  when  all  the  laser  ions  are  in  the  ground  state;  and 
D=1  when  they  are  all  in  the  excited  state. 

As  for  ae(As)  and  aa(As)  generally  we  can  measure  them  by  experiment  data  as  well  as 
McCumber  relation. 
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(3-5) 


^a(^)  =  ^e(^)exp 


h(v-£) 
k,T  ; 


where  he  represents  the  net  thermodynamical  free  energy  AF  required  to  move  one  rare  earth  ion 
from  the  ground  state  to  an  excited  state,  while  maintaining  the  lattice  temperature  T  constant. 
With  peak  value  of  aa(A,s=1.53pm)  =7e-25  m^  and  ae(7.s=1.53pm)  =  0.92aa(7.s=1.53pm).  the  Er^"^ 
density  is  assumed  to  be  p=le25  m'  . 

As  for  the  pump  power  and  signal  power  along  the  fiber  we  can  use  the  following  coupled 
differential  equation  to  describe  [Saleh,  1990  ]. 


dPpizJ)  ^  . 


dz 

dPsiz,t) 

dz 

dz 


=  +p;t^ )  +  a^p;- 
=  P,r^(a,,N,-a^M-(^sPs 

=  ±PX XeX  -  (^saN, )  ±  2a,,N,T^hu,Au  + 


(3-6) 


The  cross  sections  are  a,  the  loss  term  a  represent  internal  loss  of  the  amplifier,  T  are  signal-to- 
core  overlap  and  pump-to-core  overlap. 

Equations  (3-1)  through  (3-6)  can  be  called  numerical  model  of  EDFA  developed  by  E. 
Densurvire  et  al. 

There  is  another  popular  model  developed  by  Saleh  et  al.  [Giles,  1991]  The  simplified  version 
should  be  [Calvez,  2001]. 


er = er  {q--  q-  )-aA 

or = S  g;  expj  ^  *  {&■  -  Q""  )-a,Lj 

Where  photon  flux  Qk=Pk/hvk. 
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Fig.  3-1.  Energy  level  diagram  corresponding  to  a  basic  three-level 
laser  system,  where  the  laser  transition  occurs  between  levels  1 
(ground)  and  2  (metastable).  The  symbols  R,  W  and  A  correspond  to 
pumping  rates,  stimulated  emission  rates  and  spontaneous  decay  rates 
between  related  levels;  superscripts  R  and  NR  refer  to  radiative  and 
nonradiative  emission,  respectively. 


3.2.Simulation  results  and  discussion 

We  used  equations  (3-4),  (3-5)  to  plot  the  gain  coefficient  vs.  wavelength  at  different  D  in  Fig. 
3-2.  It  can  be  seen  from  the  figure  that  for  D  =  -1  all  ions  are  in  the  ground  state  and  the  medium 
is  absorbing  at  all  signal  wavelengths,  as  the  gain  is  negative.  As  the  relative  inversion  increases, 
however,  a  spectral  region  near  the  long  wavelength  side  of  the  transition  is  characterized  by  a 
positive  gain  coefficient.  For  these  wavelengths,  the  medium  is  amplifying,  while  for  the  rest  of 
the  spectrum,  the  medium  is  still  absorbing.  As  D  increases  toward  the  maximum  value, 
representing  complete  medium  inversion,  the  region  of  positive  gain  coefficient  widens  to  spread 
eventually  over  the  whole  spectral  range. 

Using  equations  (3-7)  and  (3-8)  we  simulated  gain  vs.  fiber  length  and  pump  power  that  are 
shown  in  Figure  3-3  and  3-4  respectively. 
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Fig.  3-2.  Signal  gain  coefficient  as  a  function  of 
wavelength  for  different  relative  medium  inversion  D 


Fig.  3-3.  Simulation  results  of  gain  vs.  amplifier  length  under 
different  pump  powers 
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Fig.  3-4.  Simulation  result  of  gain  vs.  pump  power  under  different  amplifier  lengths 
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CHAPTER  IV 


MODELING  AND  SIMUALTION  OF  RAMAN  AMPLIFIER 


4.1.  Model  description  and  simulation  results 

4.1. 1.  Forward  pumping 

The  interaction  between  pumping  source  and  signal  wave  can  be  described  by  following  coupled 
equations  [Agrawal,  2001  ]: 


dz 

£± 

dz 


(4-1) 

(4-2) 


Where  Ip  is  the  pump  intensity,  Is  is  the  signal  intensity,  gR(X)  is  the  Raman  gain  coefficient 
(material  dependent),  and  ap  are  the  losses  at  the  signal  and  pump  wavelengths.  In  general, 
gR(A,)  depends  on  composition  of  the  fiber  core  and  can  vary  significantly  using  different 
dopants. 

Decompose  the  intensities  into  axial  and  radial  components 

I  =  r^{r,e)*r\z)  (4-3) 

So  the  optical  power  can  be  expressed  as: 


P{z)  =  2;7'|  I(r,z)rdr  =<  I{r,z)  > 


Integrate  differential  equations  over  the  transverse  dimensions: 


d^ 

dz 

dz 


-a  P  p 

P  3  aR  P  s 

-«A+^P,P. 

^eff 


(4-4) 


(4-5) 

(4-6) 


where 
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(4-7) 


<!-<>>. <lf  > 

^  Trad  Trad 


^  j  raa  j  r 


> 


Assume  pump  depletion  due  to  the  stimulated  proeess  is  negligible  we  have 


dP, 

— ^  =  -a  P 
dz 


so 


-  =  -(^sPs  + 


(4-8) 

(4-9) 


We  integrated  Pg  over  z 


P,(Z)  =  P,(0)exp 


^/-,(0)V -a.L 


(4-10) 


where  =  ^—(1  -e  )  is  the  effective  length. 


a. 


I  directly  solved  the  coupled  ordinary  differential  equations  numerically.  The  calculation  results 
are  shown  in  Fig  4-1. 


Fig  4-1.  Simulation  results  of  different  input  pump  (left)  and  signal  (right)  power  evolution 
in  forward  pumping  DRA 
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4.1.2.Backward  pumping 


The  differential  equations  for  a  backward  traveling  wave  of  the  signal  is  given  by  [Mochizuki, 
1985  ] 


where 


^^^  =  [7Sp{z)-a\p^{z)  +  Y-hv^ -S^iz) 


-  L)] 


(4-11) 


(4-12) 


where  Ps/p  is  the  signal/pump  power.  Ss/p=  Ps/p/S  is  the  effective  power,  y  is  the  gain  constant. 

If  neglect  y-hv^-  Sp{z)  term  we  can  solve  the  function  analytically  and  the  solution  is  given  by 


Ps (^)  =  Ps (0) exp]  ^  •  exp(-apZ)  •  [expia^z)  -  l]-a^Z  \ 


S-a 


(4-13) 


Alternatively  we  can  solve  it  numerically.  Fig  4-2  shown  the  simulation  results. 


Fig  4-2.  simulation  results  of  different 
pumping  input  of  backward  pumping  DRA 
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4.1.3.Bidirectional  pumping 

The  models  mentioned  before  are  all  uncompleted  since  many  factors  that  must  be  considered 
in  DRA  design  were  neglected.  They  include  double  Rayleigh  scattering  (DRS),  spontaneous 
Raman  scattering.  The  most  complete  model  for  RA  design  is  given  by  [Kidorf,  1999  ].  The 
following  effects  were  included  along  with  the  capability  of  modeling  any  number  of  signals 
and  pumps  [Kidorf,  1999  ]: 

■  Spontaneous  Raman  emission  and  its  temperature  dependence; 

■  Rayleigh  scattering  including  multiple  directions; 

■  Stimulated  Raman  scattering; 

■  Amplification  of  the  above; 

■  Arbitrary  interaction  between  an  unlimited  number  of  pumps  and  signals  from  either 
direction:  Pump-pump,  pump-signal,  and  signal-signal  interactions; 

■  High-order  Stokes  generation; 

■  Combinations  of  multiple  fiber  types,  filters,  isolators,  etc.,  within  a  span. 

The  effects  that  were  excluded  are:  anti-Stokes  generation,  polarization  effects,  time 
dependence,  and  nonlinear  index  effects.  The  computer  model  is  given  by  following  equations 


=  +a{v)P^-  (z,  V)  +  r{v)P  +  (z,  V)  +  X  •  [P"  $•)  +  P"  d-P" 


K  A 

g>v  ^eff^eff 


h{g-v)IKT 


+  hv^A^-^\p\z,g)  +  P\z,g)]\\  +  - 

L 

+  Z  ^)1-  P" 


1 


-1 


■dv} 


K  A 

g>v  ^eff^eff 

Aeff 


1  + 


1 


Mv-g)/KT 


dg] 


(4-14) 


In  this  equation  (-I-)  is  refer  to  forward  propagation  direction  and  (-)  is  refer  to  backward 
propagation  direction.  Pf  (z,  v)  and  Pb(z,  v)  are  forward  propagation  power  backward 
propagation  power.  a(v)  is  the  optical  attenuation  and  y(v)  is  the  Rayleigh  scattering 
coefficient  at  frequency  v,  h  is  the  Plank’s  constant,  k  is  the  Boltzmann’s  constant,  and  T  is  the 
temperature  of  the  fiber. 

The  simplest  case  is  the  DRA  compose  of  one  forward  pumping  source,  one  backward 
pumping  source  and  one  signal  wave.  The  simulation  result  is  shown  in  Fig  4-3 
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Fig  4-3.  Simulation  result  of  a  bi-directional  pumping  DRA 
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